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We calculate the first dynamical evolutions of merging black hole-neutron star binaries that con- 
struct the combined black hole-neutron star spacetime in a general relativistic framework. We treat 
the metric in the conformal flatness approximation, and assume that the black hole mass is sufB- 
ciently large compared to that of the neutron star so that the black hole remains flxed in space. 
Using a spheroidal spectral methods solver, we solve the resulting field equations for a neutron 
star orbiting a Schwarzschild black hole. The matter is evolved using a relativistic, Lagrangian, 
smoothed particle hydrodynamics (SPH) treatment. We take as our initial data recent quasiequi- 
librium models for synchronized neutron star polytropes generated as solutions of the conformal 
thin-sandwich (CTS) decomposition of the Einstein field equations. We are able to construct from 
these models relaxed SPH configurations whose profiles show good agreement with CTS solutions. 
Our adiabatic evolution calculations for neutron stars with low-compactness show that mass trans- 
fer, when it begins while the neutron star orbit is still outside the innermost stable circular orbit, 
is more unstable than is typically predicted by analytical formalisms. This dynamical mass loss is 
found to be the driving force in determining the subsequent evolution of the binary orbit and the 
neutron star, which typically disrupts completely within a few orbital periods. The majority of the 
mass transferred onto the black hole is accreted promptly; a significant fraction (~ 30%) of the mass 
is shed outward as well, some of which will become gravitationally unbound and ejected completely 
from the system. The remaining portion forms an accretion disk around the black hole, and could 
provide the energy source for short-duration gamma ray bursts. 
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I. INTRODUCTION 



The infall of compact objects into black holes (BHs) is 
of considerable interest in many branches of astrophysics. 
In particular, many of the arguments that can be made 
about coalescing neutron star-neutron star (NSNS) bi- 
naries also apply to coalescing black hole-neutron star 
(BHNS) binaries. Both are strong candidates for the cen- 
tral engines of short-duration gamma ray bursts (GRBs), 
since the merger timescale following tidal disruption is 
comparable to the GRB duration and the gravitational 
binding energies provide the characteristic energy scales 
inferred by observers It is possible that any ejected 

matter may contribute significantly to the r-process el- 
emental abundance of the universe 3, ^ ^] . Addition- 
ally, they are expected to be among the most important 
sources of gravitational waves (GWs) that can be de- 
tected by both terrestrial laser interferometers such as 
LIGO 6], VIRGO 0, GEO jl], and TAMA % as well 
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as the proposed space-based interferometer LISA 

The key difference between the sources that can be 
observed with LIGO (and comparable detectors) and 
LISA is the characteristic frequency of the GW emis- 
sion: LISA'S characteristic frequency range falls within 
10^** — 10~^ Hz, whereas LIGO operates between 10 — 
500 Hz. Because of this, LIGO is most sensitive to the 
mergers of stellar-mass BHs, whereas LISA will observe 
more massive merging systems that involve either inter- 
mediate mass BHs (IMBHs), Mbh = 10^ - IQ-^Mq, or 
supermassive BHs (SMBHs), Mbh > IO^Mq. The for- 
mation history leading to these encounters is likely to 
involve completely different processes. 

Compact binaries with stellar-mass BHs are likely to 
be formed through typical stellar binary evolution, at 
rates that depend on parameters such as the binary mass 
ratio distribution, common-envelope efficiency, and the 
physics of supernova kicks, all of which remain some- 
what uncertain (see 0, lo. IT^ and references therein 
for a thorough review). The mass distribution of BHs 
in such systems is poorly constrained, as none have been 
observed to date, but may vary widely, spanning a range 
2Mq < Mbh < 25Mq For sufficiently tight bi- 

naries, merger will occur within a Hubble time. In these 
cases, the dissipative effects of gravitational radiation will 
cause the orbit to circularize as the binary separation 
shrinks, so that the eccentricity of the orbit is expected 
to be almost zero by the time the binary enters the LIGO 
band. Whether or not the compact object is tidally dis- 



2 



rupted by its BH companion, as well as where this would 
occur in the latter case with respect to the Innermost 
Stable Circular Orbit (ISCO), depends on both the com- 
paction of the compact object and the mass ratio (see 
Section ITT|l . 

This simple picture does not apply to compact ob- 
jects orbiting BHs with considerably higher mass. Both 
IMBHs and SMBHs are expected to reside within stel- 
lar clusters, whose dynamics will be determined by both 
stellar-BH and stellar-stellar gravitational encounters 
(scattering). Some stars will typically be scattered, either 
strongly or weakly, into the "loss cone", i.e., the volume 
of phase space encompassing orbits with sufficiently small 
periastrons that the star will be tidally disrupted before 
being kicked into another orbit by future encounters (see 
[l5| for a review of the original derivations, and 0, ^3 
for more recent work). As a result, most objects that 
enter the loss cone do so at very high eccentricity, with 
periastron distances of 5 — 50 Schwarzschild radii jlMlT^ . 
In many cases, these systems will approach the BH with 
eccentricities e ^ 0.1 |2(l |. 

GW detections from coalescence with higher mass BHs 
may yield very little information about the physics of NS 
matter, for the case of a NS falling into an IMBH, or any 
compact object (BH, NS, or white dwarf) falling into an 
SMBH with M > IO^Mq. These objects should plunge 
through the ISCO of the BH intact, since the tidal- 
disruption radius lies within the ISCO, and will likely be 
swallowed whole by the BH. For the opposite case, appli- 
cable to white dwarfs (WD) falling into IMBHs (and NSs 
into stellar-mass BHs), tidal disruption will occur outside 
the ISCO, a process we describe in detail in SecUll 

For the vast majority of its lifetime, a stellar- mass com- 
pact object binary will inspiral very slowly, such that 
it can be described by a point-mass, post-Newtonian 
(FN) treatment. FN formalisms for the adiabatic in- 
spiral epoch are now completely determined up to 3.5PN 
order |2ll |. and include lowest-order spin-orbit and spin- 
spin terms |2^. Once finite-size and tidal effects 
become important at close separation, it becomes nec- 
essary to solve the fully nonlinear Einstein field equa- 
tions. Quasi-equilibrium binary configurations in cir- 
cular orbits have been calculated in GR for NSNS 
m llpl m m El 113, BHNS (see 113, hereafter 
BSS; [SjI, hereafter TBFS, and references therein), and 
BHBH binaries ([H H Hi Hiljfor a thorough review 
of the topic and references, see |33,|33). Details of the 
transition from slow inspiral to rapid plunge, and devi- 
ations from the point-mass energy versus frequency re- 
lation found in quasi-equilibrium sequences, may yield 
important information about the physical parameters of 
the NS equation of state (EOS; see, e.g., |Mlllli3)- K 
has been suggested that a combination of 10 — 50 
broadband and narrowband observations of NSNS merg- 
ers might be able to constrain the NS radius to within 
a few percent. We will show below that BHNS mergers 
may be just as interesting, but it is likely that the inter- 
pretation of physical features in the GW signal will be 



significantly more complicated, since differences between 
stable and unstable modes of mass transfer may lead to 
radically different scenarios. 

Eventually, for those systems in which the tidal limit 
is reached before the ISCO, mass transfer onto the BH 
will begin. This process is fundamentally dynamic in na- 
ture, and can only be modeled accurately by relativistic, 
three-dimensional hydrodynamic calculations. Attempts 
have been made to model these systems analytically, but 
as we will show below, the conclusions rely on a number 
of unphysical assumptions. The earliest work describing 
mass transfer in detail for compact object binary merg- 
ers assumed that mass transfer in NSNS binaries 
would conserve both mass and orbital angular momen- 
tum, and that both NSs would remain on a quasi-circular 
orbit in corotation during the evolution; a similar set of 
assumptions was used to describe BHNS binaries as a 
possible source of gamma-ray bursts |43j . A more com- 
plex treatment developed in [44] drops the assumption of 
circularity, since it is not seen to hold in numerical cal- 
culations (e.g., 3). Still, their model for the evolution 
of BHNS systems undergoing mass transfer depends on 
a number of ad hoc assumptions that need to be tested 
by dynamical calculations in order to be proven valid. 

Beyond uncertainties about the form of the late- 
inspiral GW signal produced by a BHNS merger, there 
remains the question of the event rate, which remains 
uncertain given the complete lack of detection of such 
systems to date. Still, it is possible to estimate the likely 
merger rate using population synthesis models, which can 
be calibrated against the observed galactic NSNS binary 
population and supernova rates. Recent estimates pre- 
dict an advanced LIGO annual detection rate of any- 
where from a few mergers [T^ up to potentially several 
hundred p^ . 

Should a BHNS binary merger be observed, it might 
reveal a great deal about the physics of matter at nu- 
clear densities. In particular, the onset of mass transfer 
would yield a clear indication about the NS radius, and, 
as we will explain in detail below, the stability of the 
mass transfer would yield important information as to 
the nuclear EOS. Whereas for NSNS binaries the char- 
acteristic frequencies of GW emission during the merger 
and formation of a remnant (either a hypermassive NS or 
BH) will typically occur at frequencies outside the peak 
sensitivity of even an advanced LIGO detector, the same 
is not true for BHNS binaries. Since the frequency at 
the onset of instability scales roughly inversely with the 
total binary mass, we expect stellar- mass BHNS merg- 
ers to occur at characteristic frequencies at which LIGO 
will be most sensitive, ~ 100 — 500 Hz. If the GW sig- 
nal from a merger was observed to be coincident with 
a short-duration gamma-ray burst, we could potentially 
determine their distance, luminosity, and characteristic 
beaming angle ^ detailed theoretical understand- 

ing of these systems is now more urgent than ever, in 
ligh t of the recent localizations of short GRB afterglows 
El EllilliUHl, the first ever for these systems 
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(many long-duration GRBs have been localized, but are 
believed to be the result of collapsing stars, not merging 
compact binaries). 

Unfortunately, the current state-of-the-art for hydro- 
dynamic calculations of BHNS inspiral and merger is far 
behind that for NSNS mergers. Calculations of the latter 
have been performed using a variety of Newtonian, PN, 
and relativistic gravitational formalisms (see, |33,|53l ^^r 
thorough reviews, and 53], hereafter FGR, for a more 
recent summary) . Many calculations have now been per- 
formed in either the conformal flatne ss ( CF) approxima- 
tion to general relativity (GR) [s^, 13, or in full GR 
HE Hi, iHl, liil. These GR calculations now include 
sophisticated treatments of the NS EOS and physically 
appropriate initial spins (H^; NSs are expected to be 
nearly irrotational in the inertial frame, since the viscous 
timescale is much longer than the inspiral timescale, see 
[Eii3 and Sec. HD below). 

The key difficulty that must be overcome to perform 
simulations of relativistic BHNS mergers is the same one 
that arises in the study of BHBH binaries; the presence 
of a spacetime singularity inside the black hole. To avoid 
encountering the singularity in a numerical simulation, 
the BH interior is excised from the computational grid 
in most current applications. This is justified by the fact 
that no information can propagate from the BH interior 
to the exterior, so the exterior can be evolved indepen- 
dently of the interior. While progress has been reported, 
especially very recently |6l|, black hole evolution cal- 
culations have been plagued by numerical instabilities. 
In some ways, BHNS mergers are even more difficult to 
evolve consistently, since both the singular behavior of 
the BH as well as the hydrodynamic nature of the NS 
must be confronted. Whereas the BHBH problem in- 
volves a pure vacuum solution of the Einstein field equa- 
tions, the NS must always be evolved in such a way that 
the relativistic fluid is treated properly. 

As a result of these difficulties, all hydrodynamic calcu- 
lations performed to date of stellar-mass BHNS mergers 
have used Newtonian or quasi-Newtonian gravitational 
treatments [H H HI Ulil HI . Needless to say, bi- 
naries containing a BH can be evolved accurately only by 
using relativistic hydrodynamics in a relativistic space- 
time. We emphasize here that this applies both to the 
tidal field created by the BH, as well as the self-gravity 
of the NS. Previous Newtonian calculations have in some 
cases 0, |6^ used an approximate black hole potential, 
suggested in [b^I, that creates an ISCO at 6Mbh, but 
no single static potential can generate the full set of 
relativistic forces experienced by matter in the strong- 
field regime. Calculations employing a fixed background 
BH metric have typically been performed for stars un- 
dergoing a tidal interaction with a massive BH, rather 
than a stellar-mass BH, with relativistic dynamical terms 
but a Newtonian treatment of the self-gravity. The sec- 
ondary, in fact, is often assumed to be a white dwarf or 
main-sequence star. These models include SPH treat- 
ments without self-gravity ^^l, and both PPM ^] and 



spectral method [TOj treatments with Newtonian self- 
gravity. More recently, SPH techniques have been de- 
vised that evolve the NS matter in the background metric 
of a stellar-mass BH, using Newtonian-order correction 
to model the NS self-gravity, for both SPH [tI ^ and 
characteristic gravity |73l |. This approach is appropriate 
for describing main-sequence stars or white dwarfs. How- 
ever, since the tidal disruption is a result of a competition 
between the black hole tidal force and stellar self-gravity, 
this approach is not sufficient to describe BHNS binaries 
accurately. Modeling tidal disruption in BHNS binaries 
requires a relativistic treatment of both the black hole 
and the neutron star. 

Here, we will make use of the CF approximation to GR, 
introduced by Isenberg [t^I and Wilson and collaborators 
|75|. The CF approximation amounts to assuming that 
the spatial metric remains conformally flat, so that the 
gravitational fields can be found by solving the constraint 
equations of GR, decomposed in the conformal thin- 
sandwich (CTS) decomposition 76], alone. The GTS 
formalism has been used in numerous applications to con- 
struct initial data descri bing bo th NSNS and BHNS bina- 
ries in quasiequilibrium [Mill 111 il 111 IM IM El 113 ■ 
For these initial data the choice of a conformal back- 
ground metric is completely consistent with Einstein's 
initial value (constraint) field equations, although dif- 
ferent choices may describe the astrophysical situation 
at hand more or less accurately. The situation is dif- 
ferent for dynamical simulations in the CF approxima- 
tion (e.g. mm III 111 III]), since the assumption that 
the spatial metric remains conformally flat is no longer 
strictly consistent with Einstein's fleld evolution equa- 
tions. For many applications, however, CF provides an 
excellent approximation. For spherically symmetric con- 
flgurations, as an example, the CF approach is exact, and 
for many other applications the error has been shown to 
be in the order of at most a few percent (see, e.g., 0|). 
It is particularly useful for exploring dynamical behav- 
ior, e.g., collapse or tidal break-up, which occurs on dy- 
namical timescales and is unaltered by secular effects like 
gravitational radiation-reaction. 

In this paper we present the dynamical extension of 
BSS, who calculated the flrst relativistic, quasiequilib- 
rium BHNS sequences as solutions of the CTS decompo- 
sition of the Einstein field equations. Modifying their 
code to treat the metric for the Schwarzschild BH in 
isotropic (CF) coordinates, rather than the Kerr-Schild 
coordinates reported in BSS, we take their corotating 
quasiequilibrium configurations as initial data. As in 
BSS, we assume an extreme mass ratio, Mbh 3> -A-^ns, 
which allows us to hold the BH position fixed and re- 
strict the computational grid to a neighborhood of the 
NS, thereby avoiding complications arising in the BH in- 
terior. We also assume a polytropic equation of state for 
the neutron star, as well as synchronous rotation. The 
resulting dynamical calculations are the first of their kind 
to solve the CTS field equations for the spacetime around 
the NS self-consistently by treating both the NS and BH 
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relativistically. They allow us to study details of the 
dynamical mass-transfer process, particularly its stabil- 
ity. The CF approximation holds a stable equilibrium 
configuration constructed in the CTS formalism in strict 
dynamical equilibrium. Our calculation is a prototype 
of more detailed general relativistic calculations we hope 
to provide in the future that will involve irrotational NS 
models with more realistic EOSs and compactions, arbi- 
trary mass ratios, and a fully self-consistent treatment of 
the spatial metric. 

In the CF approximation, gravitational radiation re- 
action must be added in by hand in order to drive the 
system toward merger. While it is the secular energy 
losses to gravitational radiation that initially drive the 
binary system toward the point of tidal disruption, they 
play a much reduced role in the dynamics thereafter. In- 
deed, while secular forces determine the path the binary 
takes prior to merger, the merger itself is a fundamentally 
dynamical process, as we discuss in great detail below. 

Our work is organized as follows. In Sec. ^ we discuss 
the important physical scales that define our problem, 
and present a detailed treatment of the traditional pic- 
ture for determining the stability of mass transfer. We 
then discuss the limitations of this model, and explain 
why it may not be applicable for BHNS mergers. In 
Sec. mil we describe our numerical methods, including 
the details of both our implementation of the CF field 
equations as well as our use of smoothed particle hydro- 
dynamics (SPH) techniques to evolve the fluid configu- 
ration. In Sec. II VI we compare our relaxed initial data to 
previous quasiequilibrium models, and find that we can 
construct configurations that satisfy the field equations 
to high accuracy while reproducing previous results. In 
Sec. we present our simulations of merging binaries 
for different models of the NS polytropic EOS. Finally, 
in Sec. IVII we discuss our results in the context of GW 
astrophysics, and describe our plans for further calcula- 
tions. 



II. PHYSICAL OVERVIEW 

The evolution of binaries containing NSs is a fully rel- 
ativistic problem, since lowest-order PN approximations 
break down in the strong gravitational fields present dur- 
ing late stages of the merger. However, we can use in- 
formation from Newtonian and quasi-Newtonian calcu- 
lations to estimate the various timescales and physical 
regimes we expect to encounter. Thus, we first classify 
the relevant physical scales we expect to encounter in our 
study of BHNS binary evolution, and then generalize the 
standard model for stable, binary mass transfer to rela- 
tivistic stars. In doing so, we will explain why this model, 
variants of which have been used previously to describe 
the evolution of compact binaries, is unlikely to apply to 
BHNS mergers. 

Our simple mass-transfer model does have physical rel- 
evance, as it can apply to the case of a WD inspiraling 



in a nearly circular fashion toward an IMBH in a glob- 
ular cluster. Such a star will begin transferring mass 
long before reaching the ISCO. However, since WDs are 
typically kicked into highly eccentric orbits prior to in- 
teractions with the BH, the orbit may not have time to 
circularize fully before the onset of mass transfer. In such 
cases, the binary evolution will be more complicated than 
the scenario we consider here: it has been studied before 
by several groups (see |iall3i references therein). 

A. Units, Timescales and Characteristic Lengths 

The four most important timescales characterizing the 
problem at hand are the NS dynamical timescale tjj, 
the viscous timescale t^, the orbital timescale T, and 
the GW radiation-reaction timescale tew- Through- 
out this paper, we set G = c = 1. The BH and NS 
masses can be written in terms of the initial mass ra- 
tio q as A/bh = q^^Mj^s or equivalently Mns — qMbu, 
and the NS radius i?NS = C~^Mns — gC~^MBH, where 
C = A/ns/^ns is the compactness parameter. 

The dynamical timescale of the NS is given by 



We wish to compare this with the orbital and radiation- 
reaction timescales at the radius where Roche-lobe over- 
flow will begin. To estimate this radius, we will use the 
approximate form proposed in |8lj . 

/ \ 1/3 

Rr = 0.46a f ) ' 

which gives the (volume-averaged) Roche-lobe radius as 
a function of the mass ratio and binary separation a, for a 
binary treated as a pair of point masses in the corotating 
frame. This definition differs from the original definition 
of the Roche lobe, which was defined for incompressible 
matter (point masses are in effect infinitely compress- 
ible), but the physical scalings are the same; the coeffi- 
cient becomes 0.41 for incompressible matter instead of 
0.46. The Roche-lobe radius is equal to the NS radius at 
a separation 

= 2.17g2/3(i + g)i/3c-iMBH, (3) 
at which point the (Keplerian) orbital period is 

r^2^yj = 20.1(1 + g)i/2C-i-5A/NS 

^ 20.1(7(1 + g)i/2c-i-5AfBH - 14iD(4) 

where Mt = Mns + Mqh is the total binary mass, and 
the last relation holds for g <C 1, or equivalently, Mt — 
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-^BH- We note that in this hmit, the orbital period is a 
fixed multiple of the NS dynamical timescale, regardless 
of the properties of the NS. 

For a point-mass binary on a circular orbit, lowest or- 
der radiation reaction predicts that the binary will inspi- 
ral, losing energy and angular momentum, on a charac- 
teristic timescale icw given by 

_a _E _ J _ a* 

where r is the coalescence time, i.e., the remaining time 
until the point-mass binary would reach a = 0. At the 
Roche-lobe separation, assuming q <^ 1, we have 

^GW = 1.73g2/3c-%/Ns 

= l.TSg^/^C^^'MBH - 1.23g2/3c-5/2t^. (6) 

In general, the radiation-reaction timescale will be at 
least an order of magnitude longer than the dynamical 
timescale for any binary which begins mass transfer out- 
side the ISCO radius. The orbital period T, however, 
may become similar to the coalescence time ^gWi indi- 
cating that the infall becomes quite rapid. 

In Fig. ^ we show the regions in parameter space 
for which the critical separation for Roche-lobe overflow 
lies within or outside of the ISCO, for a wide variety of 
compact object-BH binaries. Dashed vertical lines corre- 
spond to the approximate compactness of either a WD or 
a NS, whereas dotted horizontal lines show the approxi- 
mate mass ratios to be expected for a IOA/q stellar-mass 
BH, an IMBH, or a SMBH. Systems above the critical 
curve reach the tidal limit before the ISCO, and are likely 
to transfer mass onto the BH. For those sufficiently be- 
low the curve, we expect that the compact object will 
pass through the ISCO intact, and plunge onto the BH 
relatively intact. 

We note however, that this simple picture may very 
well be altered by a number of more complicated effects. 
Recently, Miller j82i] has argued that even if systems are 
expected to reach the mass-shedding limit prior to cross- 
ing the ISCO, in many cases they will have already begun 
to plunge. Indeed, since the binary energy as a function 
of separation fiattens out significantly near the ISCO, 
tcwy which is already nearly of order T, will system- 
atically underestimate the infall velocity (a similar ar- 
gument was used in [Zp] to argue that the GW energy 
spectrum produced in NSNS mergers declines dramati- 
cally near the ISCO). Because of this, the ISCO may sys- 
tematically underestimate the binary separation at which 
prompt merger becomes inevitable. 

On the other hand, describing the "plunge" of an ex- 
tended object like a NS may provide a misleading picture 
of the dynamical merger in cases where the tidal dis- 
ruption occurs inside the ISCO but outside the horizon. 
While it is certain that some matter, likely a significant 
fraction of the NS mass, will plunge inward directly onto 
the BH, this may liberate a great deal of angular mo- 
mentum into the outer parts of the NS [s^. As a result. 
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FIG. 1: The critical mass ratio q — M,/Mbh as a function 
of the secondary's compaction C — M* / R-t , for which mass 
transfer begins at the ISCO, taken here as a = 6Mbh (solid 
line). For systems above the curve, mass transfer begins while 
the orbit is stable; for those below, the secondary may plunge 
into the BH before being tidally disrupted. Dashed vertical 
curves show characteristic compactions of a WD or NS; dotted 
horizontal curves show typical mass ratios for a IMq compact 
object orbiting a lOAf© stellar-mass BH, a lO^M© IMBH, or a 
10'' Mo SMBH. Dot-dashed curves show where f3, the ratio of 
the light crossing time the viscous timescale of the NS, equals 
unity, and where avis, the nondimensional turbulent viscosity 
assumes its largest reasonable value for the turbulent viscosity 
[See Eqs. Q and @]. Only configurations to the left of these 
curves can synchronize prior to merger. 



some fraction of the mass may survive the plunge, at least 
initially, in the form of a "mini-NS" , which will escape 
outside the ISCO on an elliptical orbit. Needless to say, 
only dynamical calculations will clarify the role played 
by these competing effects. 

The final timescale we must consider is the viscous 
damping timescale tvis, which we expect to play a crucial 
role in determining the fate of the binary once mass trans- 
fer commences. In the limit that the viscous timescale 
is extremely short (high viscosity), we expect two im- 
portant phenomena to occur. First, tidal dissipation can 
synchronize the binary so that the secondary is corotating 
upon the onset of mass transfer. From [s^, we note that 
a binary will be synchronized by the time mass transfer 
begins only if 

/3=:^ > 60(l + g)5/V'/'C3, (7) 

tvis 

Wis < i:{l + q)-'/'q'^'C-^M^s 
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< ^(l + 9)''/V/'C-4AfBH, (8) 

where /3 is the ratio of the Ught crossing time of the sec- 
ondary to its viscous timescale, as defined by jH^- If 
we follow 01 and assume that turbulent viscosity is the 
primary damping mechanism, we can define avis, the 
turbulent viscosity parameter, so that 

avis = 7^- (9) 

tVis 

We see that synchronization will occur if 

aY,,>m{l + qf'\-^/^C''\ (10) 

On Fig.^ we show curves marking the critical mass ratio- 
compactness dependence for (3—1, which we define as 
the "causal limit", as well as for avis = 0.1, which is the 
maximum plausible value for turbulent viscosity in phys- 
ical systems of interest . Configurations to the left 
of the curve can synchronize before merger; this includes 
essentially all mergers where the secondary is either an 
MS star or a WD. NS mergers, on the other hand, will 
be irrotational in general, especially when the primary is 
a BH, since the required viscosity to synchronize the NS 
increases as the primary mass increases [GflJ . Viscos- 
ity should also play a role after mass transfer starts, as 
we will discuss in detail below. 



B. The stability of mass transfer 

Once the secondary fills its Roche lobe, it will begin 
to transfer mass onto its companion. Such a process can 
be either stable or unstable, depending on its response 
to mass loss. If the volume of the Roche lobe shrinks 
faster than (or expands slower than) the stellar radius, 
the process is unstable, and the star will typically be 
disrupted violently. On the other hand, if a small amount 
of mass loss causes the star to shrink back within the 
Roche lobe, it is possible for the mass loss to temporarily 
cease, or at the very least settle down to a much smaller 
equilibrium level, whose value can be determined based 
on the assumptions made about conservation of mass and 
angular momentum, as we discuss below. 

We first note that models of stable mass transfer typi- 
cally assume that the binary orbit remains quasicircular, 
which in turn is only possible if the viscous timescale is 
short relative to the orbital and GW timescales. Main- 
taining a circular orbit requires that the orbital energy 
evolve according to a fixed relation in terms of the or- 
bital angular momentum and the mass of the secondary, 
but there is no reason to assume that such a relation 
should hold a priori. Indeed, it is viscous dissipation that 
drives the orbit toward circularity, by converting excess 
orbital energy into other forms. Furthermore, when the 
viscous timescale is long, the mass-transfer rate can grow 
extremely rapidly, since the inner Lagrange point trav- 
els into the secondary at roughly v^Rns /o-R: unbinding 



progressively denser material from the NS. If this leads 
to an unstable runaway, it is the mass loss that drives 
the orbital evolution, and we expect to find the develop- 
ment of an orbital eccentricity. This violates the typical 
assumptions made in conservative mass-transfer models, 
which assume that mass loss is steady, and slow enough 
that the orbit can remain circular as mass is lost. 

The early attempt to follow the mass transfer process 
in detail for NSNS binaries was provided by 42], who 
modeled the heavier NS as a point mass, and the lighter 
secondary using an EOS that yields a nearly flat mass- 
radius relation down to Mns ^ 0.3Mq, below which 
the NS begins to expand rapidly with further decreas- 
ing mass. Rather than assume conservative mass trans- 
fer, they parameterized the possible loss of both mass 
and angular momentum from the system, finding that 
the former has very little effect on their results. In their 
model, mass transfer leads to a widening of the binary 
orbit, under the condition that the NS radius must equal 
the radius of its Roche lobe. As mentioned above, this 
will only hold for systems in which iyis ^ ^gw- Over 
time, the mass loss rate and GW luminosity decrease 
rapidly from their large initial values at closest approach 
(as does the rate of neutrino production as the NS matter 
decompresses during the transfer), until eventually the 
low-mass NS begins to expand rapidly [s^ and unstable 
mass transfer begins. 

Many of these ideas were revisited for a discussion of 
BHNS binaries in , in light of the optical identification 
of GRB counterparts at cosmological distances. Assum- 
ing a Newtonian rt = 1.5 polytropic EOS and fully con- 
servative mass transfer, they find that the initial mass- 
transfer rate between a 1.4Mq NS and a BH with mass 
Mbh = 3 — 5Mq will occur at a rate of ~ IOOMqS"^ 
for approximately 1 ms before decaying away according 
to the approximate power-law relation A/ns oc t~^^^^^, 
which corresponds to M^s{t) oc t^^^^^. As in they 
assume that the process will terminate when the NS 
reaches a critical minimum mass and begins to expand 
unstably. 

The most recent treatment of BHNS coalescence makes 
a completely different set of assumptions about the dy- 
namics during mass transfer. Based on the Newtonian 
BHNS numerical calculations of Rosswog, Speith, and 
Wynn 0, Davies, Levan, and King |4J| assume that the 
rapid timescale for mass transfer will violate the assump- 
tion of circular orbits, which underlies the typical conser- 
vative, quasiequilibrium mass-transfer formulation. In- 
stead, they make the following assumptions: 

1. Mass transfer occurs during a timescale corre- 
sponding to half an orbit. 

2. During this time period, the NS, treated as a uni- 
form density sphere, will lose mass from a shell 
whose depth is a distance equivalent to the infall 
rate from the beginning of the mass-transfer rate 
multiplied by half an orbital period. 
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3. Half of the angular momentum lost to the trans- 
ferred mass will return to the NS, placing it on an 
eccentric orbit that will typically not lead to over- 
flow during the next periastron passage. 



This model does reproduce well the extremely high mass 
loss rates initially seen duringNewtonian numerical cal- 
culations of BHNS mergers but the assumptions 
adopted are somewhat ad hoc. In particular, transferring 
angular momentum back to the NS without adjusting its 
mass causes a discontinuous evolution of the binary orbit. 
In some cases, the NS will find itself on an orbit whose 
periastron is outside the mass-shedding limit, leading to 
a period of stable evolution until GW dissipation forces 
the orbit to decay inward again back to the onset of mass 
transfer. In contrast, we find below that mass transfer 
can be quenched temporarily, but from this point on the 
NS follows an elliptical trajectory that will take it back 
within the mass-shedding limit prior to the next perias- 
tron passage. 

In Appendix we derive a semi- analytic formulation 
for conservative mass transfer onto a BH, modeling sec- 
ondaries either by a Newtonian or a relativistic poly- 
trope. We recover the scaling relations found in 42, 43], 
and generalize them for arbitrary polytropic indices. Al- 
though these relations are unlikely to hold for merging 
BHNS systems, as shown in Fig.^ the Newtonian results 
can be applied to merging WDs as well as main sequence 
stars undergoing mass transfer. We note that there are 
semi- analytic formalisms for describing non- conservative 
mass transfer as well (see, e.g. jS2| for a formalism in- 
volving mass transfer from a main sequence star onto a 
companion), and that these have been useful in describ- 
ing WDBH mergers Ull, but that the actual NS tidal 
disruption process is sufficiently dynamic that essentially 
all analytic treatments break down. 

The theory of accretion disk dynamics presents sev- 
eral interesting connections to that of merging binaries, 
since questions about the stability of mass transfer ap- 
pear as well (see, e.g., js^ and references therein). One 
key difference between the models is the typical radial 
angular momentum distribution; parameterizing the tan- 
gential velocity profile as vt {r) cx r" , irrotational NS have 
a nearly flat velocity profile, a ^ 0, and corotating NS 
a flat angular velocity profile, a = 1, both larger than 
the Keplerian value aK = —0.5. Moreover, NS differ 
greatly from disks because of their infall velocity when 
they pass through the ISCO, and their strong self-gravity. 
While angular momentum distributions with larger val- 
ues of a help to stabilize mass transfer in disks, stronger 
self-gravity destabilizes mass transfer |^. Thus, it is 
hard to generalize across the classes, although we note 
that irrotational NS should, if anything, be more prone 
to unstable mass transfer, as the NS loses more angular 
momentum per unit mass lost from its inner edge. 



III. NUMERICAL TECHNIQUES 

To compute the dynamical evolution of a BHNS bi- 
nary, we fix the position of the BH and assume that the 
surrounding spacetime metric takes the form appropri- 
ate to a nonspinning Schwarzschild BH. The approxima- 
tion of a fixed BH position is correct in the limit that 
Mbk ^ A^NS- Here, we will study binaries with mass 
ratios q = Mns/A^bh = 0.1, which is presumably within 
the range of values for which the approximation of an 
extreme mass ratio is valid. 

To calculate gravitational forces and evolve the fluid 
configuration, we will work within the CF formalism, 
which we explain in more detail in Section IIII Bl below. 
We assume that the spatial metric is remains conformally 
flat, so that it can be written in the form 

ds^ ^ {-a^ + PkP'')dt^ + 2/3,dx'dt + i;^S,jdx'dx^ , (11) 

where a and (3^ are the lapse function and shift vector, 
respectively. Under this assumption we only need to solve 
the 3-1-1 constraint equations for ip, a and to determine 
the metric. 

Our initial configuration places the NS in a corotating 
initial configuration. Irrotational configurations, which 
are more realistic astrophysically, will be treated in a 
later publication. We model the NSs as relativistic poly- 
tropes, and assume adiabatic evolution, which we de- 
scribe in detail in Sec. IIII Al 

The code we use both to relax and evolve BHNS bina- 
ries is similar to that introduced in FGR for evolving 
NSNS binaries. We solve the five linked non- linear field 
equations of the CF formalism, Eqs. HH), ((23, and 
below, using the LORENE libraries, publicly available at 
http : //lorene . obspm . f r These Poisson-like equations 
are solved using spectral methods, decomposing the fields 
and their sources in a set of radially distinct domains into 
radial and angular expansions. Dynamical evolution is 
treated through SPH discretization. Many aspects of the 
code were discussed in detail in FGR, so we concentrate 
instead on the changes and new features introduced to 
evolve BHNS binaries. 

Roughly speaking, we have made three significant 
changes to the code to admit the presence of a BH in 
the binary. First, the asymptotic Schwarzschild BH con- 
tribution to the spacetime metric is held fixed, allow- 
ing us to solve the field equations describing the self- 
gravity of the NS in a fully consistent way. Second, as 
discussed below, we solve Poisson-like elliptic equations 
for Tp and (aV'), as in BSS and elsewhere, rather than for 
1/ = \na and (3 = ln(a-0^), as in FGR and related treat- 
ments (TBFS denotes the latter quantity "ct"). Third, 
we restrict the spatial domain of our spectral methods 
field solver to a finite radius centered on the NS, as was 
done in BSS and TBFS, which allows us to avoid prob- 
lems near the BH. Indeed, our computational domain 
is chosen so as not to overlap the event horizon at any 
time. As a result, we do not make use of the asymptotic 
boundary conditions typically used by LORENE-based 
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codes, which can be extended to spatial infinity through 
the proper coordinate transformations |91,J . The use of a 
restricted spatial domain has been introduced before, in 
the context of domains with ingoing and outgoing GWs 
|9^ . but with a set of BC's that arc not appropriate 
to the (elliptic) problem at hand. Instead, as we de- 
scribe below, we have introduced a multipole expansion 
BC, used here and in TBFS, which should be more ac- 
curate than the lowest-order power-law falloff conditions 
traditionally used in grid-based calculations. Below, we 
first summarize the relevant equations that comprise rela- 
tivistic hydrodynamics fSec. lIII X|l and the CF formalism 
fSec. lIIIBl . introduce the "split" equations which factor 
out the BH contributions to the spacetime in Sec. IIII CI 
describe our new approach for introducing a multipole 
BC in Sec. IIII Dl and finally describe how this affects 
the evaluation of various quantities in the SPH evolution 
equations Sec. IIII El 



A. Relativistic Hydrodynamics 

We assume that the matter can be described as a per- 
fect fluid so that the stress-energy tensor takes the form 

T^" = Po(l + + (12) 

Po 

where po, e, P, and denote the rest mass density, 
specific internal energy, pressure, and 4-velocity, respec- 
tively. We will describe the NS by a relativistic polytropic 
EOS that evolves adiabatically with index F. Hence, the 
pressure obeys the relation 



P = (F - l)poe, 



and initially satisfies 



(13) 



(14) 



where k is a constant. As discussed in BSS, we can 
scale away dimensional units by setting k = 1 (see their 
Sec. IIIc). 

The Lagrangian continuity equation (FGR, 54]) can be 
written as 



+ p^div"" = 



(15) 



where we define the conserved density 

p^ = au^ip^PQ ^ jri->P^po, (16) 
and the coordinate velocity 



(17) 



and introduce the Lorentz factor for the fluid 7„ = au*^. 
Lagrangian time derivatives are related to Eulerian par- 
tial time derivatives through the familiar relation d/dt = 



d/dt + v^di. To determine the Lorentz factor, we solve 
the normalization condition for the 4-vclocity, 



7„ = [au') = I + -J4- = 1 + - 
implicitly. 



Tkp 



r-i 



(18) 



The Euler equation can be written 



dUt "V'®;, p , On , ~ a flj , 2fca(7n - 1) a , 

—— = OiP - ahu^dia + UjOifr -\ ■ o^yj, 

dt -fnW 

(19) 

where the speciflc momentum is defined by 

Ui = hui, (20) 
and the specific enthalpy h by 

h=l + re. (21) 
Finally, the energy equation takes the form 



— + e^d,v' = 0, 
dt 



(22) 



where e* = jntp^ipo^^^^Y^^ ■ For an adiabatic evolution 
without shock heating, the energy equation is satisfied 
automatically by adopting Eq. (|14|l . 

To account for shocks, we included an artificial viscos- 
ity prescription composed of both linear and quadratic 
terms (the relativistic analogue of the form introduced 
in [g^l, similar to that found in Q)- We found no evi- 
dence for significant shocks within the body of the NS, as 
only the matter in the mass transfer stream directed to- 
ward the BH showed signs of significant heating very near 
the BH. Using the value of k = P/p^ = (r - l)e/Po"^ 
as a measure, a quantity that remains constant during 
an adiabatic evolution, we found variation of no more 
than 5% within the body of the NS. This is hardly a 
surprise, as there is no physical mechanism such as a 
collision to cause significant shocking within the bulk of 
the NS. Shock heating will be important for understand- 
ing the evolution of the initially low-density accretion 
stream that falls toward the BH, especially near the event 
horizon. In this region, the heating can be substantial, 
but it seems not to introduce significant feedback on the 
NS remnant. Given these results, we replace the energy 
equation, Eq. H22|) . with its adiabatic solution, Eq. I|14|) . 
throughout the calculations described here. In future cal- 
culations, where shocks may be more important, we will 
restore the full evolution of the energy equation with an 
artificial viscosity prescription and allow for shocks every- 
where. This will be especially important for irrotational 
NS calculations, since the matter transferring through 
the inner Lagrange point has significantly greater angu- 
lar momentum than in the irrotational case, and a great 
deal of it will likely forming a disk rather than accreting 
promptly. 
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B. The CF formalism 

In the CF formalism Olz^ we assume that the spatial 
metric is not only conformally flat initially, but that it 
remains conformally flat. In particular, for the 3-metric 
we approximate dtlij = so that in rectangular coor- 
dinates 7y = 5ij at all times. Strictly speaking, this 
is inconsistent with Einstein's evolution equations, but 
is often a very good approximation, particularly on dy- 
namical timescales when secular motion due to radiation- 
reaction is not important. Under this approximation the 
evolution equation for the spatial metric yields a relation 
between the extrinsic curvature and the shift, 



Kin = j- 



(23) 



Inserting this expression into the momentum constraint 
yields an equation for the shift 

S/^pi + ld\dj(3^) = Wna^P'^iE + P)W 

+2ai;^K'^\/j{ln[a/t//]) = 5^,(24) 

where is the flat space Laplacian. The Hamiltonian 
constraint is an equation for the conformal factor tp 

VV = -27rV'^£; - li^'^K.jK'^ = S^. (25) 
8 

To derive an equation for the lapse a, the remaining un- 
determined function in the metric, Eq. we choose 
maximal slicing K = ^"^^Kij = at all times, which 
implies dtK = 0. This choice can be combined with the 
evolution equation for the extrinsic curvature, which then 
yields 

V2(aV') = 2-Ka^^{E + 25) + \a^^K,jK'^ = S^^. (26) 

8 

In the above equations the matter sources E, S and [/' 
are projections of the stress-energy tensor T^^ and can 
be expressed as 



E = Poh^l-P, 



5 3P 



72 - 1 

In 
In 



[E + P), 



(27) 
(28) 

(29) 



In practice, it is easier to decompose the three coupled 
equations for the shift, Eq. into four decoupled Pois- 
son equations. To do so, we follow ^Q^, 0| and define 



/3^=4i?,--[9,(x + i?fea:")] 



and solve the set 



(30) 

(31) 
(32) 



These Poisson-like equations, found in |9J| and else- 
where, are exactly equivalent to those found in FGR for 

= Ina and /3 = ln(aV'2), and share the same asymp- 
totic fall-off behavior, but have radically different proper- 
ties near the horizon of the BH, where the lapse function 
goes to zero. This causes divergences in the values of i/ 
and /3, whereas a and remain finite and easy to deal 
with in a numerical treatment. Our chosen variables also 
exhibit a slightly different behavior when we split them 
into additive pieces contributed largely by the NS and 
BH, i.e., the contributions from the NS and BH to tp and 
atp are additive, whereas the logarithmic dependence of 
the 'V — /3" set means that the two contributions are 
combined multiplicatively. 

As several different sets of notation have now been in- 
troduced into the literature to define equivalent quanti- 
ties in the CF formalism, we present alternate notations 
used in a selection of other works in Tabled 



C. BHNS binaries 

The CF approximation is exact for spherically sym- 
metric configurations, reproducing the TOV equation for 
fluid configurations as well as the Schwarzschild solution 
for a stationary, non-spinning black hole. In isotropic 
coordinates, such a solution is given by 



2_ /l-AfBH/2r 



ds^ = - 



1 + AfBH/2r 



1 



BH 



2r 



Sijdx^dx-' . 



(33) 

From this metric we identify the BH lapse and conformal 
factors as 



aBH = 



BH 



1 ~ MBH/2r 
l + MBH/2r' 
^ ^ Mbh 
2r 



(34) 
(35) 



The BH contribution to the shift (/3i)BHi vanishes in 
isotropic coordinates (unlike in the Kerr-Schild coordi- 
nates used by BSS). 

To convert this line element to the more familiar 
Schwarzschild (areal) coordinates, with 



ds^ 



1 



2Af, 



BH 



r 



dt"- 



1 



2A/, 



BH 



r 



df' 



one makes the coordinate transformation 
r = I IH z — I r. 



2r 



1 r 



f- - Mbh + - 2AfBH) 



-r^dn^, 
(36) 

(37) 
(38) 



Note that this implies that the Schwarzschild radius and 
ISCO radius take the values 



2 A/bh 

6 AfBH 



0.5 A^BH, 
4.949 A/bh- 



(39) 
(40) 
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TABLE L A comparison of our notation for various relativistic quantities to that found in a selection of previous works using 
the CF formalism: |2^ . l5Ml5^l75ll94|l . For those cases where no unique terminology was defined, we give the simplest equivalent 
algebraic form. 



Quantity 


Here FGR [53] Gourgoulhon [27] Oechslin [§4] Wilson [JS] Shibata [94] 


Lapse 


a 


N 


N 


a 


a 


a 


Shift 




-Ni 


-Ni 




Pi 




Conformal Factor 




y/A 


y/A 




<!> 




Rest Density 




p. 




P* 




P* 


Lorentz Factor 




7n 






W 


au" 


Velocity 






NU^ + N^ 


v' 






Specific Momentum 


Ui 


Ui 


Wi 


Ui 




Ui 


Enthalpy 


h 


h 


h 


w 


h 


1 + F6 



At asymptotically large distances, f — r + Mbh- 

As discussed at length in [23|, it is useful to "split" 
the field equations when dealing with binaries, so that 
the terms on the RHS of each equation are concentrated 
on either component of the binary. Here, the method 
is slightly different. Since the BH solution is an exact 
solution of the vacuum field equations, it can be sub- 
tracted out of the full metric field equations to yield the 
CF solution for the largely-NS contribution to the fields. 
Defining N = aTp, and accordingly, A^bh = ckbhV'bh — 
1 — A/bh/2?', we split the fields such that 

V' = 'IpBH + lpNS, (41) 

a = QBH + aNS, (42) 
N = aip ^ Abh + Ans (43) 
^ Ans = aNsV'NS + "bhV'ns + ^nsV-'bh- 

The NS piece of the field equations, Eqs. and 
can be expressed as 

VVns = -27r(7/)BH + V'Ns)^£^ 

'li^BH + ^bm)HK^J)m{K'')m, (44) 

+ ^(A^BH + AfNs)(^BH + ^Ns)*(^*j)Ns(if'^)NS-(45) 

The BH contributes to the shift vector, Eq. H24|l only 
through the lapse function and conformal factor, since 
the black hole contribution to the shift vanishes in 
isotropic coordinates. 



D. Multipole Boundary Conditions 

The LORENE-based field solver we use decomposes 
the angular dependence of all scalar, vector, and tensor 
quantities into spherical harmonics (the radial decom- 
position into Chebyshev polynomials is described in de- 
tail in |9lj). For configurations in which the outermost 
boundary extends to spatial infinity, the outer boundary 
condition can be set exactly to zero for any field which 



satisfies a power-law falloff. For a BHNS binary, how- 
ever, that is not an option, since we encounter numerical 
difficulties when the computational domain overlaps the 
BH singularity. Instead, we must impose an approximate 
BC for each field on the outermost (spherical) boundary, 
which lies at a finite radius, as shown in Fig. ||2Jl. This 
outermost boundary is chosen so that it never overlaps 
the BH event horizon. 

Any Poisson-like equation V^<i> — p with compact sup- 
port has an exterior solution given by 



$(f,0,, 



OO I 

EE 

OO I 



= E E PlmXlmiO , (j)) 



(46) 



1=0 m=-l 



where we define the multipole moments of the source 
term pim [see Eq. (4.2) of jo^]- This solution estab- 
lished the boundary conditions for our outermost com- 
putational domain, and can be matched to the interior 
solutions to yield the field solution everywhere in space. 
Here, the source terms of the Poisson-like Eqs. (|24ll . H25|l . 
and H26|l are not compact, but instead satisfy rather steep 
power-law falloffs, allowing us to use the same formalism 
while introducing only small errors. Noting that the mat- 
ter configurations are equatorially symmetric, we only 
sum over multipoles with the same equatorial symmetry 
as the particular field (i.e., / -I- m even for tp, a, /3^, 
and x'j I + ^ odd for /3^). Rather than evaluate the real 
field source integrals against the complex spherical har- 
monics Y/„i , we evaluate both the multipole moments and 
the resulting expansions against the real and imaginary 
parts of the spherical harmonics Y/m with m > (noting 
that 1/0 are purely real and that Y;,-™ = (-1)™1'/J„)- Fi- 
nally, we truncate the expansion at a predetermined value 
I — Imax, where throughout this paper we use Imax — 4, 
or hexadecapole order. Thus, we assume pim = for all 
terms with I > Imax when we define the BC's for our 
field equations. This is done for two reasons. First, the 
multipole coefficients fall off steeply at high /, so that the 
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higher order multipole make ever smaller contributions 
to the field at large separation. Second, including higher 
order multipoles can lead to purely numerical instabilities 
in the field solvers for a finite set of Chebyshev polyno- 
mials, since the rapid oscillations with respect to angle 
can lead to large gradients in derivative-based quantities. 

We find that that a multipole treatment can lead 
to significantly higher accuracy for our boundary solu- 
tion, at the cost of some computational efficiency. To 
avoid numerical instabilities arising from quickly growing 
higher-order multipoles, we employ underrelaxation dur- 
ing each iteration, updating each field such that ipnew = 
(1 — X)tpoid + X4'new, where ipoid is the field value from 
the previous iteration, and ipnew is the new solution found 
from solving the elliptic equation. We find good stability 
and efficiency by setting A = 0.5 initially, and increasing 
the value to A = 0.5 — 0.05 log(A/3^) with each itera- 
tion, where A/J" is the maximum relative change in the 
y-component of the shift vector from iteration. The iter- 
ation loop terminates when < 10~^, at which point 
A ~ 0.95, representing very weak underrelaxation. 

Our multipole BC's allow us to calculate the forces 
on particles that fall outside the computational domain 
directly, since the multipole expansion for the metric is 
valid throughout space. Indeed, for such particles, we 
calculate the BH contribution to the lapse and confor- 
mal factor from Eqs. H34|) and (|35|) . the NS contribution 
from the multipole expansions given by Eq. (|46|l . and the 
gradients of the NS contribution from 



9$ 



^21 + 1 



E 

l.rn 



Pin 



„2l+l 



„2l+3 



-(47) 



noting that dr/dx^ = x^/r. 

This works directly for the lapse and conformal factor, 
but the shift vector is slightly more complicated. Recall 
that we have solved elliptic equations not for but for 
Bi and ^ defined in Eq. (|30|l . and thus only know 
the multipole decomposition of the latter quantities. In 
terms of these, the gradient of the shift is given by 



(48) 



where we also evaluate terms of the form 
9^$ ^ fd^djir'Yi^) 



dx^dx^ 



l,m 

21 A 



^21 + 1 



1 



^21+3 yx'd,{r'Yi„,) + x^d,{r'Yi„^) 

+S,jr^Yi„i) 
{2l + l){2l + 3)x'x^r%rn 

r2l+5 



(49) 



Since the lapse goes to zero at the horizon, particles ap- 
proaching it become frozen in proper time, and cannot 
penetrate within. 

The approach we use has several advantages over the 
leading order power-law falloff BC's typically used in BSS 
and other grid-based field calculations (e.g., [E^)- First, 
we lose less information about the source terms by ex- 
tending to higher order multipoles ( 57] include dipole 
order falloff terms for the lapse and conformal factor in 
full GR, while [H^l include quadrupole-order terms for 
these in CF gravity). Moreover, we avoid a problem asso- 
ciated with symmetries present in our quasi-equilibrium 
initial conditions which are broken during the dynami- 
cal evolution. In particular, as we show in Appendix IbI 
our quasiequilibrium configurations can be shown to have 
a vanishing monopole contribution to and vanishing 
monopole and dipole contributions to [3^. Once the bi- 
nary becomes tidally disrupted, however, we expect the 
monopole contribution to and the dipole contribu- 
tion to to grow in magnitude {P^ may never have a 
monopole contribution, since equatorial symmetry holds 
for dynamical configurations as well) . While these terms 
are growing, we are faced with a situation where the 
leading-order term may very well not be the largest mag- 
nitude multipole contribution on the boundary. Defining 
a global power-law falloff index to fit the boundary condi- 
tion, as in previous treatments, is impossible even when 
the two lowest-order moments arc known, since the in- 
dex varies with angle. Instead, our multipole summation 
handles this situation naturally, calculating all low-order 
moments accurately. 



E. SPH Discretization, Computational Domains, 
and Timestepping 

Many of the methods used to perform an SPH dis- 
cretization of the CF hydrodynamic and field equations 
are discussed in FGR, so here we summarize briefly the 
fundamental aspects of the SPH treatment and the new 
features present in our BHNS code. The neighbor find- 
ing algorithms used in our code are based on routines 
from StarCrash, a publicly available, extensively docu- 
mented Newtonian SPH code, which can be found online 
at www . astro . northwestern . edu/theory/StarCrash. 

Motivated by the form of the Lagrangian continuity 
equation H15|l . we define the mass nia of each particle, 
fixed in time, in terms of the conserved density p*, such 
that 



{p*)a = rribWab 



b 



(50) 



where Wab is the C^, piecewise, 'W4" smoothing kernel 
function for a pair of particles introduced by [93, and 
used in FGR and elsewhere. For each particle, we de- 
fine a smoothing length ha, and compute all sums over 
particles that lie within a sphere of radius 2ha surround- 
ing each particle (we calculate all SPH quantities using 
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a "gather-scatter" technique, as described in FGR and 
the StarCrash documentation). Smoothing lengths are 
updated using underrelaxation in order to maintain a 
roughly constant number of neighbors for each particle, 
set at the beginning of each run. Each particle is ad- 
vanced through space with a velocity — dx^ /dt, which 
we evaluate with a second-order accurate leapfrog evo- 
lution scheme, calculating forces from the Euler equa- 
tion H19(l at the half-timestep. Since the calculation is 
adiabatic, the energy equation, Eq. H22I) is automatically 
satisfied when we use the adiabatic EOS, Eq. H14|l . A 
typical timestep in our evolution scheme, started with 
particle velocities evaluated half a timestep in advance 
of the particle positions, involves a number of compu- 
tational elements. First, we advance all particles a full 
timestep, and re-evaluate the particle neighbor lists and 
the SPH expressions for the density of each. We then use 
the SPH form for the density at each particle position to 
define the computational domains used by the Lorene 
field solver, shown schematically in Fig. [3 To do so, we 
calculate the position of the NS center of mass from all 
particles having a density {p<t)a > Pcrit, where Pcrit is a 
critical value chosen to encompass the vast majority of 
particles at the beginning of a run. Next, we calculate the 
surface of the innermost computational domain i?i (0, 0) , 
as the smallest triaxial ellipsoid, centered on the NS cen- 
ter of mass, that contains all particles that lie at greater 
radii from the BH than the NS center of mass, treat- 
ing the particles as spheres of radius 2ha. This is very 
similar to the technique described in FGR, except that 
there we included all particles that passed the density 
cut, regardless of which side of the NS they fell within. 
Here, however, the dynamics of the mass transfer are dif- 
ferent. In equal-mass NSNS binaries, the NS only begin 
to disrupt at very close separations, never deviating par- 
ticularly far from an ellipsoidal configuration up to the 
point of merger. Here, mass transfer is initially one-sided 
toward the BH, and the outer half of the NS remains 
virtually intact while the inner half becomes deformed 
by the tidal gravity of the BH. We find that our field 
solver performs best if we define our elliptical domain 
based on the profile from the outer half of the NS, as it 
can handle without difficulty field sources located outside 
the innermost domain, but produces numerical errors if 
the density field of the NS drops to zero within the in- 
nermost domain, (see the second panel of Fig. [21). The 
two outer domains, which have the topology of spherical 
shells, are defined initially such that their outer bound- 
aries are spheres at radii equal to twice and three times 
that of the maximum extent of the innermost shell, i.e. 
R2{0,(t)) = 2 X max(i?i);i?3(6',0) = 3 x max[i?i]. Over 
time, we hold the outermost boundary fixed at this ra- 
dius, and adjust that of the second domain to be the geo- 
metric mean of the outer radius and the maximum value 
from the inner domain, i.e., R2 — 0.5 * {R3 + max[i?i]) 
(compare the two panels of Fig. [21). 

Once the computational domains are defined, we use 
the techniques of (oH to define a set of "collocation 
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FIG. 2: A schematic representation of the spectral methods 
computational domains used to solve the NS components of 
the field equations, Eqs. I25ll - 124|l . Initially, a triaxial ellipsoid 
with surface Ri{6,(f>) and origin at the NS center of mass is 
fitted around all SPH particles for which (p*)i > Pcrit (top 
panel) . Two annular "shell- like" domains with spherical outer 
boundaries are also laid down, with radii R2 and R3, twice and 
three times the maximum value of R\. During the evolution 
(bottom panel), we use the same procedure to fit Ri(6,(j}), 
keep the value of R3 fixed, and calculate R2 as the mean 
of R3 and the maximum of Ri . Some particles first leave the 
innermost domain, and then the entire computational volume, 
particularly those accreted by the black hole, shown as a circle 
of radius 0.5 Mbh centered at the origin. 



points" at which we compute the local SPH expression 
for , u, and P* = Kp^ , noting the latter remains equiv- 
alent to Eq. H13|) . for adiabatic evolution and polytropic 
initial data. From these, we calculate all other hydrody- 
namic quantities using the Lorene library routines, and 
solve the field equations iteratively. After every iteration 
of the field solver, all hydro quantities are updated to 
reflect the new fields. 

Once a convergent solution is found, we must export 
back all relevant matter and field terms from the spectral 
decomposition to the particle positions. For particles in 
the innermost domain, we evaluate most hydrodynamical 
terms directly from the spectral decomposition. Thus, 
denoting by "SB" those terms evaluated in the spec- 
tral basis and "SPH" those quantities defined only on 
a particle-by-particle basis, we calculate the Euler equa- 
tion as 

f duA g \diP' 

\ / m L J SPH 
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ahu'^dia 



2ha{^l - 1) 



dill) 



SB 



SB i"jJsPH- 



(51) 



This approach works in the outermost domains for ex- 
trinsic quantities hke that go to zero smoothly at the 
surface of the NS matter, but fails for intrinsic quantities 
that have discontinuities there, , e.g., and 7„, since the 
Chebyshev radial decomposition cannot describe discon- 
tinuous functions. Instead, we evaluate hydrodynamic 
terms for particles in these domains on a particle-by- 
particle basis, and evaluate field quantities and deriva- 
tives through the spectral decomposition. 



dui 
~dt 



KISB 



SPH 



-[a9ia]sB[/iw°]spH 



SB 



In 



SPH 



SB NJSPH- 



(52) 



After calculating the forces for the RHS of the Eulcr 
equation, we advance the velocities from their original 
half-timestep value forward to a half-timcstep ahead of 
the positions, and then resolve the field equations to 
determine the velocity using the same approach de- 
scribed above for particles based upon their computa- 
tional domain, 



1 



[ui]sPH - [/3'']sB, 



(53) 



SB 



in the innermost domain, with vP evaluated via SPH in- 
stead for the outer ones. 



IV. EQUILIBRIUM MODELS 



SPH value for the density of each particle, and adjust 
the masses and smoothing lengths of each particle until 
each has approximately the correct number of neighbors 
as well as the correct density, to within ~ 2%. While the 
resulting configuration could serve as acceptable initial 
data, we can do better by evolving the configuration in 
the corotating frame with drag forces applied, to damp 
away spurious deviations from true equilibrium. This 
also allows us to relax to quasiequilibrium initial models 
with binary separations differing by up to ~ 20% in ei- 
ther direction using the same initial data from BSS. Of 
course, the new field solution will be different, reflect- 
ing the change in magnitude of the tidal terms, but we 
have found that after approximately 1000 timesteps of 
relaxed evolution, the overall level of spurious motion is 
equivalent. 

We used two grid-base datasets to generate our initial 
data. For configuration A, the NS is modeled as a rel- 
ativistic F = 1.5 polytrope, of compaction Mns/-Rns = 
0.042 (or equivalently, mass Mns = 0.05 orbiting a BH 
of mass Mbh = IOMns at a separation oq = ILSA/bh- 
Configuration B features a NS with the same compaction 
but a stiffer EOS, F = 2 and oq = ll.lAfBH (where M 
is the dimensionless mass defined in Sec. IIIC of BSS). 
Note that in these units the maximum compaction of an 
isolated NS is 0.214 and 0.074 for adiabatic indices F = 2 
and F = 1.5, respectively. 

To convert from the units of BSS to those used here, 
many quantities must be linearly rescaled. In particular, 
for configuration A, f = 10.5, and M-qu — 4.72. Thus all 
distances should be multiplied by a factor r/M-Qu = 2.2 
to convert from the "hatted" units of BSS to those here 
expressed in terms of the BH mass. Similarly, for con- 
figuration B, f = 1.32 and A/bh = 0.5, so the rescaling 
factor is 2.64. 



The first step in evolving BHNS binaries is the con- 
struction of accurate initial data. In our approach, this 
requires not only determining the fields and hydrody- 
namic quantities within and surrounding the NS, but also 
the construction of a relaxed SPH discretization config- 
uration describing the NS itself. 

We take as our starting point data constructed from 
the grid-based equilibrium scheme described in BSS. We 
modified the scheme of BSS to allow for a conformal back- 
ground metric corresponding to a Schwarzschild BH in 
isotropic coordinates, which can be adopted more eas- 
ily for the CF approximation used here, rather than the 
Kerr-Schild background used in BSS (see also TBFS). 
To construct an SPH particle decomposition, we first lay 
down a hexagonal close-packed lattice of SPH particles 
over the Cartesian coordinate volume where the density 
of the star is positive. Tentative particle masses are as- 
signed to be proportional to the density p*, normalized 
to match the proper NS mass. Next, we calculate the 



Both NS models are undercompact compared to the 
expected physical NS parameters. Since our method as- 
sumes an extreme mass ratio, we are limited to low com- 
pactness NS models in order to study cases where tidal 
disruption occurs outside the ISCO, as can be seen from 
Fig. ^ Thus, while our configurations do not exactly 
represent physical parameters expected to be found in 
BHNS binaries, they serve as an analogue to binaries 
containing lower-mass BHs and more compact NS that 
will have comparable tidal-disruption radii located out- 
side the ISCO. In a future work, we will treat more phys- 
ically realistic NS compactnesses, as well as NS spins, 
including cases for which the tidal-disruption radius is 
within the ISCO. 

Below we describe the technique by which we generate 
our relaxed SPH initial conditions in Sec. II and then 
show the comparison between our resulting models and 
the grid-based data in Sec. IIVBI 
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Relaxation of Initial Data 



which impHes 



When preparing a fluid configuration to be evolved us- 
ing SPH, it is generally necessary to use some form of 
relaxation first. Otherwise, numerical deviations from 
equilibrium present in the discretized initial configuration 
will drive the dynamics, leading to a variety of spurious 
effects. Relaxation is easiest to perform for configurations 
in which the matter will be stationary in some reference 
frame, such as a corotating system, since the spurious 
component of each particle's velocity can be easily iden- 
tified and damped away by a drag term in the force equa- 
tion. This statement holds equally true for Newtonian 
and relativistic formalisms, although the latter require a 
slightly more complicated numerical treatment, for rea- 
sons we discuss below, primarily due to the presence of 
velocity-dependent forces as well as a more complicated 
set of variables used to define the equations of motion. 

In order to derive the proper equations for a relaxation 
scheme in a relativistic setting, it is useful to start with 
a brief review of how the process works in Newtonian 
physics, and then generalize to the appropriate relativis- 
tic equation. In what follows, we define our coordinates 
such that the x-axis corresponds to the line connecting 
the centers of mass of the two objects, the y-direction to 
their orbital velocity, and the z-direction to the binary's 
angular velocity. 

In Newtonian physics, we have a set of inertial frame 
evolution equations 



df 

(ft 

dv 

'dt 



(54) 
(55) 



where the RHS of each are known functions used to define 
an initial condition. For the case of a corotating equilib- 
rium binary configuration, these quantities take the form 



= X r , 



(56) 
(57) 



where we use the subscript "eq" to indicate the relaxed 
value, and where fcyi is the "cylindrical" radius. To 
evolve the fluid during the relaxation, we shift to the 
frame in which the matter is stationary. Thus we define 



V 



= V — n X 



(58) 



so that V — for equilibrium configurations, and evolve 
dx/dt — V. We determine 51 as an eigenvalue from the 
condition that the binary center-of-mass separation is 
already known, by summing over all of our SPH par- 
ticles. Based on symmetry considerations, only the x- 
component of the equation yields a non-trivial result: 



(60) 



We add to the force equation a linear drag term with 
some characteristic timescale r in order to damp away 
the spurious motion in of the initial condition. The relax- 
ation timescale should generally be approximately equal 
to the dynamical time of the system. Thus, we evolve 

dV ,dv 2^ \ ^ r^l^ ^ / \ 
-TT = (-77 + fcyl) + a-drag ^ a + il Tcyl . (61) 

at at T 

As we approach equilibrium, both the term in parenthe- 
ses as well as the drag term separately approach zero. 

The relativistic case is slightly more complicated, but 
we can derive analogous relativistic expressions for all of 
our Newtonian ones. Our equations of motion now take 
the form 



dx 
du 



(62) 
(63) 



where the velocity variables are related by Eq. H17|l . v — 
u/{%Ij'^vP) ~ (3. In a corotating frame in equilibrium, we 
know that d{4''^vP) / dt = 0, and will treat the term as a 
constant. The shift vector has the time-dependence of 
the other vector quantities. 

For a corotating configuration, we have 



f2 X r , 



(64) 



and the slightly more complicated Euler equation 



dt 



= iP\'\-n^rcyi + ^x /3](65) 



The matter will be propagated on trajectories with ve- 
locity V = v — Veq, just as before. Now, however, we have 
the condition 



V — V — V, 



eq 



(3-nxr, 



(66) 



whose time derivative in equilibrium must satisfy the con- 
dition deq/{ip'^u^) — X f3 + il'^r'cyi — 0. In the particle 
decomposition, we find 



^"^'[-dT 



dV 



— j =0 = ^m,a^-^mi(a. 



eqli 



rriifl^Xi, (59) 
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which can be solved for Q. 

To the force equation, we also add a linear drag term 
with a characteristic timescale r, but the drag term must 
force u toward its equilibrium value Ueq = 4''^u^{flxr+l3), 
rather than toward zero, so that 

— — a ~ aeq = a + ip u (II Vcyi — il x p) 

. (68) 

This equation should, barring any physical instabilities, 
damp away spurious motion and produce a corotating 
equilibrium configuration. 

B. Comparison with Other Quasi-Equilibrium 
Sequences 

Our grid-based models, based on the scheme described 
in BSS but with isotropic background coordinates, were 
constructed using 48 x 48 x 24 grids, with outer bound- 
aries placed X = ±3,?/ = ±3, z — 3 for the F = 2 case, 
and X = ±2,y = ±2,z = 2 for the T = 1.5 case. SPH 
configurations were generated corresponding to these bi- 
nary separations, as well as wider and narrower binary 
separations constructed by translating the NS to the ap- 
propriate position. The number of particles used to con- 
struct the SPH configurations were Up = 103, 953 and 
Hp = 77,908 for the P = 2 and P = 1.5 EOS, respec- 
tively. For the F = 1.5 EOS configuration only particles 
of mass rrii > 10~*mi-rnax were accepted, where rni-max 
is the maximum mass of any SPH particle present. This 
mass cut is useful for eliminating an outermost layer of 
negligible total mass, which is typically blown off the sur- 
face of the NS anyway by even a tiny amount of spurious 
motion resulting from deviations from pure equilibrium 
in the initial condition. 

Each SPH configuration was relaxed for 1000 time- 
steps, which corresponds to ~ lOto, a sufficient time 
given that the initial models were rather relaxed to begin 
with. Parameters for our grid-based models and the final 
SPH configurations at the end of relaxation are listed in 
Table ini Models for F = 1.5 are labeled Al - A7, while 
models for F = 2 are labeled Bl - B4. As a check, we 
compare the value of the period determined during the 
SPH relaxation process T, to the exact relativistic Kepler 
relation for a point mass about a BH, Tn = ^/r^ /Mbk, 
and find very good agreement. Here the binary separa- 
tion is measured in ureal coordinates, whose relation to 
CF coordinates is given by H37|) . 

We note that configurations A1-A3 do not settle down 
completely during relaxation. In each case, the central 
density of the NS dropped monotonically as the configu- 
ration expanded in the x-direction, indicating that mass 
transfer would eventually begin even with drag forces ap- 
plied. 

In general, we find very good agreement between the 
grid-based initial data and our SPH configurations for 



TABLE II: Parameters for our relaxed initial models. Tn is 
the exact relativistic Keplerian period for a point mass about 
a BH, defined by Eq. 



Run 


a/MsH 




T/Mbh 


Tn/Mbh 




r = 


1.5, q = 


0.1, C = 


= 0.042 


Al 


10.438 


0.0258 


243.8 


243.8 


A2 


10.745 


0.0247 


254.0 


253.7 


A3 


11.256 


0.0232 


270.6 


270.3 


A4 


11.513 


0.0225 


279.1 


278.8 


A5 


11.767 


0.0218 


287.5 


287.3 


A6 


12.027 


0.0212 


296.5 


296.1 


A7 


12.791 


0.0194 


323.1 


322.5 




r = 


= 2,q = 


0.1, C = 


0.042 


Bl 


10.539 


0.0255 


246.4 


247.0 


B2 


10.961 


0.0242 


260.2 


260.7 


B3 


11.093 


0.0238 


264.6 


265.0 


B4 


11.648 


0.0222 


283.0 


283.3 



stable configurations. In Fig. O and 0] we show a com- 
parison between the field values and densities from our 
SPH configuration and the grid-based data along the x- 
axis, for configurations A5 and B3. In both cases the 
relevant fields agree to generally within about 1 — 2%. 
The only exception is configuration A5, where we find 
some disagreement between the two methods on the half 
of the NS facing the BH. The discrepancy is primarily 
due to the different BC's: the grid based data imposes 
a l/r power-law falloff condition on a cube whose inner 
edge is located at x = —2, whereas the multipole solution 
used for the SPH configuration is imposed on a spheri- 
cal boundary at f = 4.5. Thus, the spectral methods 
solution integrates over a much larger volume of space, 
and allows for higher-order terms in the field solution at 
the boundary, which are not insignificant at a distance 
corresponding to a few NS radii. The small disagreement 
in the density profile is in part an SPH effect: SPH typ- 
ically smooths out the density field over each particle's 
smoothing length. Since our particles are initially equally 
spaced along a lattice, this length is ^ 5x — 0.05, and 
we cannot fully resolve the sharp density peak at the NS 
center. 

We can formulate an independent check on the self- 
consistency of our initial data by checking how well they 
satisfy the integrated Euler equation, 

h 

— = constant. (69) 

This condition is typically used to generate initial data in 
grid-based calculations, e.g., Eq. (42) of BSS and Eq. (28) 
of TBFS, but appears nowhere in our relaxation scheme. 
As we show in Appendix |21 it can be derived as a con- 
sequence of a relaxed initial configuration for which the 
RHS of the Euler equation H19|l is zero. In Fig. \E\ we 
show on a particle by particle basis the value of hjvP, 
for configurations A5 and B3, plotted for clarity against 
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FIG. 3: From top to bottom, the values of p*, , ip, and 
a along the x-axis for the SPH (solid line) and grid-based 
(dashed line) data representing configuration B3, a NS with 
a r = 2 polytropic EOS and an initial binary separation 
ao/MsH = 11.1. The agreement is generally to within 1 — 2%. 



FIG. 4: The values of p*, /3^, ^j), and a along the x-axis for 
the SPH (solid line) and grid-based (dashed line) data repre- 
senting configuration A5, a NS with a F = 1.5 polytropic EOS 
and an initial binary separation ao/MsH = 11.8. Conventions 
are the same as Fig. |21 



the particle's radial coordinate position outward from the 
NS center of mass. In both cases, we find the standard 
deviation from the mean is < 0.01%, and the maximum 
discrepancy < 0.1%. 



V. EVOLUTION OF BHNS BINARIES 

From our relaxation results, it appears that the tidal 
disruption limit for the adopted choices of NS EOS would 
occur at binary separations ao/AfeH ^ 11.0, in line with 
the predictions of Eq. (0) . For models with smaller binary 
separations, we were unable to find a convergent solution 
for the configuration with a stable central density maxi- 
mum. These results can be confirmed through dynamical 
calculations in the strict CF formalism which ignore all 
energy and angular momentum losses from gravitational 
radiation-reaction. From our discussion in Sec. IIIBI we 
might expect the possibility of qualitatively different be- 
havior for NSs with the stiffer vs. softer EOS evolved 
from an initial configuration near the stability limit. For 
the softer EOS, we expect unstable mass transfer: the NS 
should disrupt completely once mass transfer begins. For 
the stiffer EOS, we might expect stable mass transfer in 
the strongly viscous regime. However, as there is no dis- 
sipative mechanism, such as viscosity, powerful enough 
to circularize the orbit after the onset of mass transfer, 
the picture is considerably more complicated. 

All models A (for F = 1.5), or all models B (for F = 2), 



describe the same physical binary system, at different 
separations representing different moments in its evolu- 
tion. Clearly, for each binary there is only one correct 
inspiral history. In our separate runs we pick up this 
history at different points (approximating the orbit as 
circular), which helps to locate the onset of tidal dis- 
ruption and analyze the system's dynamical evolution. 
Ideally, we should start an evolution calculation at some 
large separation and evolve it forward to complete coa- 
lescence, since the assumption of quasicircularity is vio- 
lated in an ever increasing fashion by inspiraling systems. 
However, even when we augment our CF equations with 
a radiation-reaction potential to drive the secular inspi- 
ral, we have neither the time nor the numerical stability 
to calculate an evolution for an indefinite time. Thus, 
our models started from different initial separations rep- 
resent a series of approximations to the true binary evolu- 
tion, which illustrate the dynamics of tidal break-up, and 
should not be taken as physically distinct evolutionary 
paths. GW energy and angular momentum losses drive 
the BHNS binary toward coalescence. However, the GW 
timescale is much slower than the dynamical timescale, 
so tew ^ ^D, we expect that radiation- reaction losses 
will cease to play an important role in the hydrodynam- 
ical evolution once phenomena associated with the dy- 
namical timescale, such as tidal break-up and mass loss, 
begin. Nevertheless, our evolution calculations started 
from outside the stability limit and including the effects 
of GW radiation-reaction yield our best models for the 
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FIG. 5: Value of the SPH expression for the integrated Euler 
equation constant hjiP as a function of the particle's radius 
from the NS center of mass for configurations B3 (top) and 
A5 (bottom), featuring a NS EOS with T = 2 and V = 1.5, 
respectively. Note that the proper value differs between the 
two cases. The standard deviation in both cases is < 0.01%, 
with maximum variation < 0.1%. 



physical evolution of the systems in question. 



A. BHNS mergers with a soft EOS: V = 1.5 

As a soft NS EOS, we choose a polytropic model with 
adiabatic index F = 1.5 (or equivalently, polytropic in- 
dex n = 2). As we see from Sec. lA 21 a NS of com- 
pactness C = 0.042 is expected to undergo unstable mass 
transfer in the high viscosity limit regardless of the bi- 
nary mass ratio (and thus in the inviscid limit as well). 
To test out how well this statement applies in the invis- 
cid limit, which applies to our calculations (see Tables 
V and VI of [o^l for numerical estimates of the viscosity 
present in a lower-resolution implementation of our cur- 
rent SPH scheme) as well as to physically reasonable NS, 
we evolve binary BHNS configurations from a number of 
initial separations. This also allows us to estimate the 
critical separation marking the onset of mass transfer, 
which according to Eq. (O should be at aji = IIMbh- 

In Fig. El we show the evolution of run A5, at an 
initial time t = 0, corresponding to the initial relaxed 
configuration, as well as tjP = 1, 2, and 2.5. The NS 
revolves clockwise around the BH, which is fixed at the 
origin. The event horizon, located at i?BH = 0.5AfBH, 
is shown as a circle. In the first plot, the NS has essen- 
tially filled its Roche lobe, and has primary axis ratios 



a^jai ~ as/oi = 0.86. We note that the figure shows 
particle locations, rather than a surface density represen- 
tation. In fact, particles near the edge of the NS have 
a density four orders of magnitude lower than in the NS 
center. After a full orbit, we see in the second panel of 
Fig. El that the NS has begun to shed a small amount of 
mass, which indicates that it is near the mass-shedding 
limit, not necessarily past it. Our initial SPH configu- 
ration is relaxed to the point where the spurious motion 
resulting from deviations from equilibrium is small, but 
not zero. Given that the dynamical timescale of the ex- 
tremely low-density outer layers of the NS is so long, 
and the SPH particle masses so small (roughly propor- 
tional to the density) , even a tiny error in the initial data 
may result in significant spurious velocities in these lay- 
ers over time. For an isolated NS, these particles will 
remain bound, but the same is not true for a NS in a bi- 
nary. Here, particles that escape the NS surface will often 
travel outside the Roche lobe and be lost from the NS. As 
a result, the NS will lose a very small amount of mass and 
angular momentum. By the third panel of Fig. El the NS 
has expanded to the point that matter is now lost through 
both the inner and outer Lagrange points. This leads to 
the formation of a stream of matter thrown out into an 
extended halo around the binary system, most of which 
remains bound to the BH. Meanwhile, there remains a 
mass stream of material accreting directly onto the BH. 
We believe that the relativistic nature of the BH gravita- 
tional potential plays an important role in the dynamics 
of this accretion process. The matter streaming through 
the Lagrange point passes sufficiently close to the BH to 
fall well within the ISCO on its first passage. As a re- 
sult, most of it accretes directly onto the BH, rather than 
forming a disk. The instability of orbits near the BH is 
likely to play an important role in suppressing the for- 
mation of an accretion disk. We note, however, that our 
assumption of an extreme mass ratio may bias the evolu- 
tion towards prompt accretion (as does the assumption 
of initial synchronization), since the BH is a fixed target, 
rather than one orbiting the binary center of mass itself. 
Finally, in the last panel, the NS is nearing a state of 
complete disruption, and will continue to do so until we 
can no longer locate a gravitationally bound object. 

A similar pattern holds for all runs we performed with 
the same soft EOS, regardless of the binary separation. 
Due to the nature of the instability, all models we cal- 
culated led to the eventual tidal disruption of the NS, 
since there is no stabilizing mechanism to suppress mass 
loss once mass transfer begins, no matter how small the 
mass-transfer rate. It does take longer for the NS to be 
disrupted in configurations placed at a greater initial bi- 
nary separation, however. In the bottom plot of Fig. 
we show the evolution of the mass loss over time for all 
of the configurations using the soft NS EOS, defined as 
the total mass that can be found outside the innermost 
computational domain at any given time. We see that 
in all cases mass loss is an unstable process, occurring 
at a rate that grows extremely rapidly until the major- 
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FIG. 6: SPH particle configurations at times t = Q, 286, 566, and 693, projected into the orbital plane, for run A5 (for which 
the orbital period is T = 287 Mbh)- These configurations correspond to the initial configuration, and after 1, 2, and 2.5 full 
orbits, respectively. The BH is shown as a circle at a radius r = 0.5Mbh. We see that once mass transfer begins, the NS begins 
to expand, forming a low-density single-armed spiral pattern. Note that particles have different masses, and that those in the 
center of the NS are in most cases significantly more massive than those originally in the outer layers. 



ity of its original mass is no longer bound to the NS. In 
the case started from the largest separation, slightly over 
half of the NS mass is accreted directly onto the BH or 
into orbits that lie within the ISCO, while half is lost 
outward to form the "spiral arm" pattern seen in Fig. 
This pattern is to be expected, as the primary response 
of the NS to losing mass is expansion. The inner half of 
the NS is pushed inward toward the BH, while the outer 
half expands outward. For a synchronized configuration, 
we know that the average specific angular momentum 
of the NS is generally greater than the specific angu- 
lar momentum near the center-of-mass of the NS, since 



j = V X r ^ flr^ is weighted more heavily by matter on 
the outside of the NS. In response to the expansion of the 
NS, we expect the orbit to tighten, which is confirmed by 
the numerical results. In Fig.^l we show the binary sep- 
aration over time for the runs with a soft NS EOS. We 
see that the onset of mass transfer leads to a rapid de- 
crease in the binary separation, prompting the explosive 
mass loss. Only in the late stages does the remnant of 
the NS begin to move back outward, but by then tidal 
disruption is inevitable. 

As the initial separation is decreased, the qualitative 
behavior of the system stays essentially the same, but 
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FIG. 7: Binary separation (top panel) and total mass lost 
from the NS (bottom panel) as a function of time for all runs 
calculated using the soft NS EOS, runs A1-A7, which have 
r = 1.5. The properties of the initial configurations, which 
differ only in their binary separation, are shown in Table HTl 
Here, particles are defined as "lost" when they lie outside the 
innermost computational domain around the NS. In all cases, 
mass loss never quenches once it begins, eventually leading to 
the complete disruption of the NS. 

more of the NS mass ends up being transferred inward 
toward the BH. In Fig. |H1 we show the amount of mass 
lost inward (top panels) and outward (bottom panels) for 
all of the runs we performed with the soft EOS. For con- 
figurations near the stability limit, we find almost twice 
as much mass falls inward toward the BH as is expelled 
outward into the spiral arm. As these calculations were 
performed without radiation-reaction effects, we expect 
that including them would tip the balance of the mass 
transfer further toward mass accretion onto the BH, since 
the binary orbit will be driven by radiation to smaller 
separations than we find here. 

Unlike the situation found in NSNS binaries, for which 
> 99% of the matter typically remains bound to the sys- 
tem, we find that a significant amount of matter is un- 
bound from the system during all of our calculations with 
the soft EOS, representing in each run between 3 — 5% of 
the original mass of the NS. This fraction is much larger 
than that typically found in relativistic calculations of 
synchronized binary NS systems [ssL Is^ . but roughly 
consistent with previous results from Newtonian BHNS 
calculations However, we note that this fraction 

is almost certainly an overestimate, perhaps greatly so: 
irrotational configurations suppress the amount of mass 
that will become unbound (see, e.g., or FGR for a 
similar argument with respect to NSNS mergers) . The 



FIG. 8; Total mass lost inward toward the BH (top panel) 
and outward away from the BH (bottom panel) as a function 
of time, for the runs shown in Fig. |7| In general, the smaller 
the initial separation, the more mass that falls inward toward 
the BH. For the runs started from a larger separation, much 
of the initial mass loss results from spurious numerical of low 
mass particles out of the Roche lobe, leading to nearly equal 
flows directed inward and outward. In contrast, the mass 
loss for the closer cases is in exactly the form expected for 
Roche-lobe overflow through the inner Lagrange point. 



total angular momentum of an irrotational configuration 
is less than that of a synchronized configuration, and the 
decrease in specific angular momentum is largest at the 
outer edge of the NS. Thus, we expect that while some 
mass may still be unbound when the NS is initially irro- 
tational, the amount may be significantly less. 



B. BHNS mergers with a stifFer EOS: F = 2 

To study a configuration that would be predicted to 
undergo stable mass transfer in the classical conserva- 
tive quasiequilibrium scenario, described in Sec. Ill Bl we 
model the NS with a stiffer, T — 2 polytropic EOS, which 
has been used in numerous studies as a first approxi- 
mation to a stiff NS EOS. We note that in Newtonian 
physics, this is the critical polytropic index for which an 
isolated NS has a radius independent of its mass. In rel- 
ativistic gravity, self-gravity effects cause the NS radius 
to increase as the NS mass decreases. 
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1. Evolution without GW radiation-reaction 

In Fig. |51 we show the binary separation (top panels) 
and NS mass fraction lost (bottom panels) as a function 
of time for runs B1-B4. We find three qualitatively differ- 
ent behaviors over time, which we will describe in turn. 
We note first that the tidal limit separation is not drasti- 
cally different for this EOS compared to the softer EOS. 
Both cases show very good agreement with the classical 
Roche limit result, as we would expect. 
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FIG. 9: Binary separation (top panel) and total mass lost 
from the NS (bottom panel) as a function of time for the 
runs calculated using the stiff NS EOS without GW radiation- 
reaction, runs B1-B4, which have F = 2. Conventions are as 
in Fig. |7| For run Bl, we see periodic mass loss, as the NS 
gets kicked into an elliptical orbit, losing mass during every 
periastron passage. Runs B2 and B3 are similar to those 
with the softer EOS, as the NS gets disrupted during the first 
mass loss phase. Run B4 is essentially stable, only showing 
the long-term effects of numerical diffusion of particles. 

Run B4 is completely stable for approximately three 
orbits. Eventually, the NS begins to lose a small amount 
of mass, which gives us an estimate for the length of time 
over which the code can reliably maintain an equilibrium 
configuration in the absence of relaxation {t/T ^ 3). We 
note from Fig.^Jthat the mass loss from the NS is almost 
evenly divided between a component moving toward the 
BH and the component directed away, indicating that 
the mass loss is a numerical artifact caused by particles 
near the edge of the NS diffusing outside the Roche lobe 
over time. 

Runs B2 and B3 are started from a binary separation 
approximately equal to the mass-shedding limit. In both 
cases, the NS makes approximately one orbit after mass 
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FIG. 10: Total mass lost inward toward the BH (top panel) 
and outward away from the BH (bottom panel) as a function 
of time, for the runs shown in Fig. ^ Conventions are as 
in Fig. |H] Approximately 50% more mass is lost inward for 
runs B2 and B3, started from near the stability limit, whereas 
for run Bl, marked by an elliptical orbit and periodic mass 
transfer, nearly 4 times much mass is lost toward the BH. 



transfer begins without any appreciable change in the bi- 
nary separation. By this point, approximately 10% of the 
NS mass is stripped away, and the binary begins to move 
outward. As we saw in the case of the softer EOS, the 
mass loss is unstable and the mass-transfer rate grows 
rapidly. In both cases, the NS is tidally disrupted before 
the mass transfer halts. By the time the NS is disrupted, 
approximately 60% of the original mass has fallen inward 
toward the BH, as we see from Fig. ^| Of the matter 
shed outward, approximately 0.075AfNS is unbound from 
the system. This result fits in with the general picture 
derived from NSNS binaries, in which mass shedding out- 
ward is more efficient for stifFcr EOS (see, e.g., 99j and 
references therein). 

The evolution of the NSs in runs B2 and B3 started 
from initial conditions differing only in their binary sep- 
arations (run B3 starts from 0.5% further out), which 
are nearly equal to the mass-shedding limit. We find 
that the closer the NS is when mass transfer begins, the 
more mass is lost inward toward the BH relative to that 
lost outward, until the NS expands to the point that it 
greatly overfills its Roche lobe and mass loss becomes 
more isotropic around the axis describing the NS veloc- 
ity. 

These observations go a great deal toward clarifying 
the evolution seen in run Bl, the only case in which the 
mass transfer was periodic, rather than continuous. Here 
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the initial configuration places the NS within the mass- 
shedding limit, as GW losses would be expected to do 
for physical BHNS systems. Indeed, runs B3a and B3b, 
described in Sec. IVB 21 which include the effects of ra- 
diation reaction, are driven to approximately this binary 
separation by GW radiation-reaction energy and angular 
momentum losses. When mass transfer begins, it occurs 
only toward the BH, which causes the NS to move out- 
ward. This expansion in the orbit, which happens with 
no outwardly directed stream to counter the outward ac- 
celeration, occurs sufficiently fast that the Roche-lobe 
expansion is enough to quench mass transfer. At this 
point, the NS is on a continuously expanding eccentric 
orbit (e ^ 0.1), whose apocenter lies outside the Roche 
limit and pericenter within it. As the NS crosses over 
the new mass-shedding limit, mass transfer begins again, 
pushing the NS onto an even wider eccentric orbit, simi- 
lar to the pattern seen in ^ for NS with a stiffer EOS. 

This scenario is ostensibly similar to that proposed by 
01) but with one crucial difference. In both cases, mass 
transfer forces the NS onto a highly elliptical orbit. Here, 
however, mass transfer occurs during every periastron 
passage, whereas in their model the NS can be kicked into 
such a widely separated orbit that GW radiation-reaction 
must drive the NS back toward the mass-shedding point. 
We believe that the assumptions that go into the latter 
model lead to this unphysical result. In |4J|, it is assumed 
that the NS loses a specified amount of mass during mass- 
transfer events but recovers half of its angular momentum 
during the next half orbital period. This can lead to a 
discontinuous evolution of the binary separation. Here, 
we see that once mass transfer ceases, the NS will follow 
a nearly unperturbed elliptical orbit, with essentially no 
change in its orbital angular momentum. Mass transfer 
must resume when the orbit crosses this same point as it 
approaches pericenter during the next passage. 



2. Mergers including GW radiation-reaction 

While evolution calculations lacking GW radiation- 
reaction terms can be useful for studying the processes 
that control the moment to moment dynamical evolu- 
tion of the system, we know that the effects of radiation- 
reaction must play an important role in the secular dy- 
namics of the system. Indeed, given the potentially un- 
stable nature of mass loss, we expect the inwardly di- 
rected component of the NS velocity to be critical. Since 
the mass transfer leads not to an instantaneous change 
in the NS velocity, but rather in its acceleration, there 
will be a time period immediately after the onset of mass 
transfer during which the NS loses mass while falling fur- 
ther inward. All the while, the inner Lagrange point will 
move further within the NS, regardless as to how its ra- 
dius adjusts on the dynamical timescale. 

To model radiation-reaction, we add a damping force 
to the material, representing the lowest-order quadrupole 
contribution to the radiation reaction potential. Thus, 



we add to the RHS of Eq. (|19|l an acceleration term of 
the form 



(70) 



where x is a quadrupole radiation-reaction potential (see 
Sec. 36.1 of (lOflj), defined here as 



(71) 



in terms of the fifth time derivative of the quadrupole 
moment, 



Q 



kl 



STF 



PABMXkXld X 



(72) 



Here padm is the quantity that can be integrated to give 
the matter contribution to the ADM mass padm = ip^E, 
which appears in the field equation for ip, Eq. (|25|l . Fol- 
lowing the argument found in Sec. Ilia of FGR, we eval- 
uate the fifth time derivative of the quadrupole moment 
in terms of the expression for the first time derivative 
in the limit that it is dominated by terms representing 
orbital motion and not changes in the density with time 
(i.e., dpADu/dt is neghgible), 



l)kl 



= STF 



PABuiXkVl + XiVk)d^X 



(73) 



where "STF" means the symmetric trace-free component 
of the tensor. We then assume that further time deriva- 
tives result purely from the orbital motion, and approxi- 
mate 



kl, 



(74) 



where the angular velocity lu is taken as the ratio of the 
angular momentum to the moment of inertia. 



{Vy)a - ya{Vx)a) 



(75) 



While the resulting radiation-reaction force will differ 
slightly from the true quadrupole expression found by 
taking the exact time derivatives, it is sufhciently accu- 
rate for our purposes here, generally within 10%. It is 
important to remember that the radiation-reaction force 
drives the binary inward on a secular timescale icwj but 
plays almost no role once effects occurring on the much 
more rapid dynamical timescale become dominant. 

To test how radiation-reaction effects change the sce- 
nario described above, we calculated two runs that in- 
cluded radiation-reaction terms. Both took as initial data 
the configuration used also for run B3, using the stiffer 
r = 2 NS EOS, placed nearly at the stabihty hmit. For 
run B3a, we added a radiation-reaction acceleration term 
described by Eq. (|7n|l : for run B3b, we doubled the mag- 
nitude of this force. The evolution of the former, run 
B3a, is shown in Fig ^2 Mass transfer occurs in a well- 
defined stream for the first two orbital periods, until the 
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FIG. 11: SPH particle configurations at times t/MsH = 0, 261, 514, and 767, projected into the orbital plane, for run B3a, 
which includes the dissipative eflects of radiation-reaction (for this configuration, T = 265Mbh)- In terms of the binary orbit, 
these correspond to the initial configuration and 1, 2, and 3 full orbits, respectively. Conventions are as in Fig. El Here, we 
see that radiation-reaction initially drives only an inwardly directed mass-transfer stream onto the BH, until somewhere after 
i/MBH = 700, when the expansion of the NS becomes unstable and tidal disruption occurs. 



expansion of the NS eventually drives the rapid tidal dis- 
ruption of the NS. 

In Fig. El we show the evolution of the binary sep- 
aration (top panels) and NS mass loss (bottom panels) 
for runs B3a and B3b. We find that the stronger the 
radiation-reaction losses, the greater the initial mass loss 
from the NS, since the first passage within the mass- 
shedding limit takes it closer to the BH. This in turn 
drives the NS back outward, quenching the mass loss tem- 
porarily until the next periastron passage, after which the 
NS tidally disrupts. From Fig. El we see that stronger 
radiation-reaction losses favors a larger amount of mass 



lost inward onto the BH (and thus less outward into a 
disk), as one would expect. Thus, we conclude that the 
inclusion of GW radiation-reaction terms have the effect 
of increasing the chance that some fraction of the origi- 
nal NS mass will remain bound after an initial phase of 
mass loss, since the NS rebounds more sharply outward 
than for cases in which radiation-reaction losses are ig- 
nored. Tidal disruption, while still seemingly inevitable, 
also occurs at a greater distance from the BH. 

In the bottom panel of Fig. El we show the gravity 
wave signal produced in run B3a, in both polarizations. 
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FIG. 12: Binary separation (top panel) and total mass lost 
from the NS (bottom panel) as a function of time for the 
runs calculated using the stiff NS EOS and dissipative GW 
radiation-reaction effects, runs B3a and B3b, which have F = 
2. Conventions are as in Fig. |7| The only difference between 
the two runs is the magnitude of the radiation-reaction effects; 
we use the quadrupole order form in run B3a, and double its 
strength for run B3b. We find that doubling the radiative 
drag force forces the binary to a smaller separation, leading to 
a larger initial burst of mass loss and a more rapid increase in 
the binary separation. As a result, the system tidally disrupts 
slower than the case shown in run B3a. 

The two components are defined by the familiar relations 

Doh+ = Qx^ - Qyy, (76) 
I?o/ix = 2Q,j,, (77) 

where Do is the distance from the observer to the bi- 
nary. The corresponding angular frequency of the GW 
signal, approximately equal to twice the orbital angular 
frequency, is shown in the top panel of the figure. Prior 
to disruption, the GW waveform takes the classic point- 
mass form, with a steadily but extremely slowly increas- 
ing amplitude and frequency (relativistic and finite-size 
affects cause minor deviations from the point-mass form; 
see ^3 for a discussion). Once the mass transfer begins, 
however, the frequency reaches a maximum and begins 
to decrease quickly, as does the amplitude. In general, if 
the NS does transfer sufficient mass to survive the initial 
infall, we expect this decrease in frequency and ampli- 
tude until the signal can no longer be reliably observed. 
Should the orbit be elliptical, as we found for run Bl, 
this will show up in the GW waveform as well. 

As matter accretes onto the BH, it will excite quasi- 
normal ringing modes that could in theory be visible in 
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FIG. 13: Total mass lost inward toward the BH (top panel) 
and outward away from the BH (bottom panel) as a function 
of time, for the runs shown in Fig. |n| Conventions are as 
in Fig. |S] Initially, in run B3b with twice the physical GW 
radiation-reaction force applied, all mass lost in run B3b is 
directed inward. This leads to a rapid increase in the binary 
separation and slows the growth of the mass-transfer rate. 
In contrast, in run B3a, mass loss is more evenly balanced 
between inwardly and outwardly directed flows, and the NS 
disrupts more quickly. 

the gravitational wave signal. Unfortunately, our numer- 
ical approach limits our ability to determine the gravita- 
tional wave signal we expect from these modes. Indeed, 
such a calculation would require a dynamical treatment 
of the spacetime very near the BH, whereas in our calcu- 
lation, the key physics for ringdown occurs in the asymp- 
totic region located outside our computational domain, 
where the BH remains stationary. To evaluate quasi- 
normal mode ringing, it will be necessary to relax the 
CF approximation, in which no gravitational radiation 
is generated, and to evolve the fields everywhere. Ad- 
ditionally, the use of a lapse function that penetrates 
the horizon will be crucial for studying the problem self- 
consistently. We plan to study these issues in future 
work. 



VI. SUMMARY AND DISCUSSION 

We have performed relativistic calculations of BHNS 
mergers in CF gravitation, in the limit that the BH is 
much more massive than the NS. These calculations mark 
the first time that a relativistic treatment has been ap- 
plied both to the self-gravity of the NS as well as the BH 
spacetime. 



24 




300 400 600 800 1000 

1/Mb„ 



FIG. 14: Gravity wave emission angular frequency (top panel) 
and waveform in both polarizations, h+ (solid curve) and hx 
(dashed curve), defined by Eqs. 17611 and 17711 . for the merger 
calculated in run B3b. Initially, we see conclusion of the stan- 
dard "chirp" signal, in which the binary separation decreases 
while the GW amplitude and frequency increases, all of which 
happens extremely gradually on the secular GW timescale. 
This lasts until the onset of mass transfer, at which point we 
encounter a much more rapid "reverse chirp" , as the GW am- 
plitude and frequency rapidly decrease while the NS is tidally 
disrupted. 



For systems studied here in which the onset of tidal 
disruption occurs outside the ISCO, previously proposed 
analytical models do not properly describe all the com- 
plex tidal phenomena that pertain to the evolution. In all 
the runs we calculated, mass transfer plays a leading role 
in determining the dynamics of the system. In general, 
mass transfer in a stream directed toward the BH causes 
the orbital separation to increase, in many cases quite 
dramatically. Mass transfer also causes the NS radius to 
expand on the dynamical timescale, with relativistic self- 
gravity effects leading to a more rapid expansion than 
is seen in Newtonian gravity for a given NS EOS. As 
a result, mass transfer is significantly more unstable in 
relativistic gravity. We conclude immediately that pre- 
vious models of mass transfer in compact object binaries 
that assume the orbit remains quasi-circular j42l l43j are 
not applicable here. Furthermore, the model put forward 
in M4| also seems to be insufficient, in that the orbital 
parameters evolve discontinuously from one orbit to the 
next (as a result of angular momentum being added to 
the NS while its mass is held fixed). We find instead 
that if some remnant of the NS survives the initial burst 
of mass loss, it can end up on an elliptical orbit that 
takes it back outside its mass-shedding separation. Dur- 



ing every successive periastron passage, however, more 
mass will be lost, eventually leading to the complete dis- 
ruption of the NS. 

As the NS expands during mass loss, it eventually loses 
mass outward as well as inward, so long as the plunge 
does not take it too far within the ISCO, leading to a 
prompt plunge onto the BH. While the majority of mat- 
ter released through the outer Lagrange point remains 
bound to the BH in the former case, a significant frac- 
tion is ejected with sufhcient velocity to become unbound 
from the system completely, approximately 5 — 7% for 
the r = 2 EOS we considered. This mass loss also limits 
the radial expansion of the orbit, and as a result we find 
in many of our calculations that mass transfer is never 
quenched once it begins. Even though the NS moves out- 
ward, it persists in configurations for which the Roche 
lobe lies within the star, and mass loss continues until 
the NS is completely disrupted. 

We plan to improve our simulations and relax several 
approximations in the near future. In particular we plan 
to adopt the astrophysically more realistic irrotational 
initial configuration of TBFS instead of the corotating 
configurations used here. Relativistic NSNS [s^,!!^ and 
Newtonian BHNS ^ calculations have shown that for 
these cases the amount of matter ejected from the binary 
system is expected to be significantly smaller, since the 
material on the outer edge of the NS has significantly 
less specific angular momentum. Thus, while we expect 
that BHNS mergers will expel more material than in the 
case of NSNSs, in which the binary components are more 
comparable in mass, we would assume that the ejected 
fractions found here are overestimates. 

Calculating mergers using irrotational NSs should also 
increase the probability that escaping matter forms an 
accretion disk around the BH, even if that disk is short- 
lived (see, e.g., 83]). We typically found in our calcu- 
lations here that most of the matter transferred toward 
the BH ends up accreting onto it directly, since the spe- 
cific angular momentum is not sufhcient to create a disk. 
Even though an irrotational NS has less total angular 
momentum than a corotating one, the matter on the in- 
ner edge has a higher specific angular momentum, and is 
more likely to orbit around the BH rather than infall di- 
rectly. In addition, prompt accretion of matter may also 
be inhibited slightly by a moving BH orbiting the binary 
center of mass, whereas in our assumption of an extreme 
mass ratio, the BH position is fixed. 

In our future work, we will test out how shock heating 
in the accretion disk affects its evolution, and determine 
if there are cases in which feedback onto the NS will affect 
its future evolution as well. By including a relativistic ar- 
tificial viscosity treatment, we will follow the thermody- 
namic evolution of the disk, as well as that of the NS and 
any outward mass loss. Of course, to investigate BHNS 
mergers fully and accurately, we will ultimately have to 
abandon the assumptions underlying the CF metric as 
well. While our description of an isolated Schwarzschild 
BH is exact, the BH lapse goes to zero just outside the 
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event horizon. This causes matter to "pile up" around 
the BH, since the proper time ceases to advance in this 
region. While this poses no difhculty for determining the 
fate of the material, which will clearly be accreted, it 
does act as a computationally challenging ridge of grow- 
ing mass concentration as material is transferred contin- 
uously toward the BH. 

These calculations adopted "undercompact" NSs be- 
cause we are interested in systems that disrupt outside 
the ISCO while restricting our attention to the case of 
extreme mass ratios. By performing dynamical evolu- 
tions using more compact configurations, we will investi- 
gate the transition to the opposite case, in which plunge 
begins while the NS is still completely bound. Based 
upon our results, we expect that that in some cases, the 
mass-transfer rate may prove sufficient to kick the core 
of the NS back out to a wider, highly elliptical orbit. 
The phase space for which this will occur, in terms of the 
NS EOS and binary mass ratio, remains poorly under- 
stood, and may not conform to simple analytic estimates, 
which have difhculty describing the unstable processes 
that characterize the merger. In it was found for a 
quasi-Newtonian potential that for a mass ratio q = 0.1, 
a NS described by a F = 2 EOS was disrupted during 
the initial passage, whereas one with a F = 3 EOS led 
to a punctuated mass-transfer scenario similar to the one 
we describe for run Bl above in Sec. IV B'2l Probing the 
assumed forms of the NS EOS that lead to complete dis- 
ruption versus survival of a remnant NS core may prove 
to be a crucial diagnostic tool for determining the true 
physical NS EOS. 

In many prior works, it has been assumed that so long 
as ajsco < o,Roche, the NS will necessarily be tidally 
disrupted, with some of its mass deposited into a disk 
around the BH. According to ^85], this may not be true. 
Instead, the plunge may start well outside of the ISCO, 
since, based on the properties of quasiequilibrium mod- 
els, the inspiral time scale may approach the orbital time 
scale already well outside of the ISCO (so that the tran- 
sition through the ISCO is nearly dynamical rather than 
adiabatic). It is unclear exactly what role the spin of 
the BH will play in this process, though it appears that 
a prompt merger is more likely for a Schwarzschild BH 
than a spinning Kerr BH j71|, where t he la t ter c ase is 
favored by binary evolution calculations |lOll Il02j | . 

Two different factors can mitigate a prompt merger, 
however, and will need to be investigated numerically in 
more detail. First, the onset of a "plunge" also marks the 
point at which we assume the quasi-equilibrium formula- 
tion to break down. From that point onward, the NS will 
no longer follow the trajectory predicted by quasiequi- 
librium results, and may therefore not plunge as fast 
as quasiequilibrium models would predict (in fact, these 
predictions provide completely unphysical overestimates 
of the infall velocity at the ISCO itself). Perhaps more 
importantly, there is no guarantee that a plunge phase 
will lead to the entire NS being swallowed by the BH 
|83l |. Angular momentum can be transferred outward 



within the NS on something approximating the dynami- 
cal timescale, and it is possible that some fraction of the 
mass, perhaps a significant fraction, will survive the ini- 
tial plunge phase. Using the techniques developed here, 
and initial configurations taken from TBFS, we will study 
how varying the NS spin and compaction affect the final 
fate of the NS. 

It will be necessary to perform relativistic merger cal- 
culations in full generality, in order to drop the assump- 
tion that the binary mass ratio is extreme. We are cur- 
rently constructing such quasiequilibrium data, which 
will then serve as initial data for dynamical simulations 
|l03j . We expect that these dynamical simulations will be 
plagued by the same difhculties encountered in dynamical 
simulations of BHBH binaries, and therefore anticipate 
that this will be a very challenging project. 

Until that point, however, a great deal can be accom- 
plished. First and foremost would be to identify the 
boundaries in phase space that separate qualitatively dif- 
ferent phenomena that can occur during a BHNS merger. 
The most obvious categories would be prompt merger, 
prompt tidal disruption leading to an accretion disk, or a 
period of periodic mass-transfer bursts, if the latter does 
occur at all. The former distinction should prove useful 
for understanding any potential X-ray/gamma-ray emis- 
sion from these systems, and aid in our understanding of 
short gamma-ray bursts like the recently observed GRB 
050509b (46, 47j', GRB 050709 M^M^ GRB 050724 50J, 
and GRB 0508I3[^. 
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APPENDIX A: MASS TRANSFER IN THE 
VISCOUS REGIME 

Below, we derive a formalism that can be used to de- 
scribe conservative quasiequilibrium mass transfer, for 
those cases where the viscosity is high, and the binary 
orbit remains circular during the mass-transfer process. 
For typical (high-mass) NSs, this regime does not apply 
([53,E3; see also Fig. ^ above), but can apply to low- 
mass NSs, WDs, and main sequence stars orbiting BHs. 
We nevertheless refer to the star as a NS below. 
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1. Newtonian polytropes 

In Newtonian physics, the mass-radius relationship is 
given as a function of the polytropic index n by 



D ji f{l — n) I (3— n) 



(Al) 



For a given value of the parameter k in P = Kpg^^^^"\ 
the familiar result is that NSs with n — 1 have a uniform 
radius independent of their mass. NSs with n < 1 have 
radii which decrease as the mass decreases, whereas those 
with n > 1 increase in size as they lose mass. We see 
that if the NS loses mass at a rate m = —M-^s (implying 
m > for the case of interest), the change in the radius 
is given by 



R 



NS 



71 — 1 m 



i?Ns 3 - n Mns 



(A2) 



Once mass transfer begins, the NS will shed mass, and 
in doing so, lose both energy and angular momentum. 
This will affect the binary orbit, in a way which depends 
on both the magnitude and fate of the angular momen- 
tum of the ejected material. 

For conservative mass transfer, the orbital angular mo- 
mentum J and the total binary mass Mt = Mns + 
Mbu = (1 -f q)M-[^s/q are conserved globally, as mass 
is transferred from the NS to the BH. We assume that 
the orbit remains circular. Since q = M^s, / [Mt — Mns), 
we see that 



q , , m 

- = -(1 + 9)77— • 
q AiNs 

From the angular momentum of a circular orbit. 



we find that 



J = MnsMbha/-^, 



[A/ns(Mt-Mns)]" 



(A3) 
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The Roche-lobe radius Rr changes during the process as 
well. Taking the logarithmic time derivative of Eq. (O 
and combining with Eq. ljA6|l . we find 
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Combining this final expression with Eq. l|A2p . we see 
that mass transfer will be unstable if i?Ns/-RNS > Rr/ Rr, 



or equivalently, 
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or in terms of the adiabatic index 
transfer is unstable if 
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In particular, the critical mass ratio for unstable mass 
transfer for some polytropic indices commonly used in 
numerical calculations are 



= 1/2 (F = 3) 
n 1 (F = 2) 
3/2 (F = 5/3) 



q > 14/15, 
q > 5/6, 
q > 2/3. 



(A12) 
(A13) 
(A14) 



As a general rule, for polytropic indices n < 3/2, mass 
transfer is stable only for systems with components of 
similar mass. For n > 3/2, the critical mass ratio drops 
quickly, down to the limiting case n = 9/4, the softest 
polytropic EOS for which stable mass transfer can ever 
occur, at which point q = 0. 



2. Relativistic polytropes 

Relativistic polytropes are not self-similar for a fixed 
value of the polytropic index; as the mass decreases, non- 
linear gravitational effects become weaker, and the star's 
scale-free density profile grows in size. It is straight- 
forward to incorporate this into our discussion of mass 
transfer. We recast the mass-radius relationship in the 
form 



-Rns 



^^'•-'ns 



7(C), 



(A15) 



where ^ sets the physical scale of the mass-radius rela- 
tion, C = M / R\s the compactness of a spherical star, and 
/(C) accounts for the relativistic corrections to the mass- 
radius relation. Since relativistic corrections effectively 
increase the strength of gravity, /(C) must be a mono- 
tonically decreasing function of compactness. As C — > 0, 
^ approaches the proper Newtonian value, — (^n, and 
/(O) — > 1. Taking a logarithmic derivative of Eq. (|A15|) 
shows us that 
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But we know 
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so we conclude that 
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3~n 3~nfdC\ fdC 
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(A18) 

From Eq. (|A8|I . we see that the critical mass ratio for 
instability becomes 



_ 9 ~ 4n 1 Cdf f Cdf 
9-3n 3-nf dC\ f dC 



(A19) 



where the second term is always negative, indicating that 
more compact configurations are more unstable against 
mass transfer. 

In Figure [T^ we show the critical compactness values 
separating stable and unstable mass transfer for relativis- 
tic polytropes, as a function of the binary mass ratio and 
the polytropic (bottom panels) and adiabatic (top pan- 
els) indices. As a general rule, as the compactness of the 
NS increases, mass transfer is more likely to be unstable. 
The Newtonian curves (C = 0) have an analytic form 
given by Eqs. (jASjl and (|A10|I . Also shown are the two 
models we evolve dynamically in this paper, both with 
C = 0.042. The case with n = 1{T = 2), shown as a 
triangle, would be expected to demonstrate stable mass 
transfer in the highly viscous regime, but as we show in 
Sec. IV Bl the true situation is nowhere near this simple 
for typical NSs. The case with n = 2{T = 1.5), shown as 
a square, is expected to lead to unstable mass transfer. 

Representing the same results against compactness, 
rather than mass ratio, shows how little parameter space 
is available for unstable mass transfer in the viscous limit. 
In Figure El we show the critical binary mass ratio for 
stable or unstable mass transfer as a function of the NS 
compactness and polytropic (bottom) or adiabatic (top) 
index. We see that for a given compactness, there is 
a rather limited range of polytropic indices which can 
produce unstable mass transfer. The heavy lines show 
the maximum possible compactness allowed for a given 
polytropic EOS; no model to the right of those curves 
can be constructed. The models that we evolve dynam- 
ically are shown as well. As noted in Sec. IV Bl many 
models expected to undergo stable mass transfer in the 
viscous limit are extremely unstable during mass transfer 
when viscosity is not a dominant driver of the evolution 
[59I l60l | , so the true parameter space for instability is ac- 
tually much larger for typical NSs than analytic models 
would otherwise predict. 



3. Time evolution 

We describe here a crude treatment of the dynamics of 
mass transfer in the presence of gravitational radiation- 
reaction losses, using a number of simplifying assump- 
tions. Most important among these will be treating the 
problem in the quasi-equilibrium regime. 
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FIG. 15: Critical compactness for determining the stability 
of mass transfer from a relativistic polytropic NS companion 
as a function of the binary mass ratio, defined in terms of the 
polytropic index n (bottom), or the adiabatic index F = 1 -|- 
1/n. Mass transfer is unstable to the right of the curve, stable 
to the left. The square and triangle represent the positions of 
the two models we are evolving dynamically. 



We will study only conservative mass transfer, in the 
highly viscous limit. The evolution of the binary sepa- 
ration can be derived from the expression for the total 
angular momentum for a circular orbit, Eq. (|A4I) . Mak- 
ing no assumptions about the evolution of the system 
angular momentum, we find 



J 
J 



Ml 



NS 



BH 



Ml 



NS 



M 



BH 



a 
2a 



M 



NS 



Setting the RHS equal to zero gives the standard expres- 
sion for conservative mass transfer, Eq. (jA6|) . For the 
case of a relativistic binary, we know angular momentum 
will not be conserved. Instead, angular momentum is ra- 
diated away from the system in gravitational waves. As 
we are assuming circular orbits, we will make use of the 
standard angular momentum loss rate. 



J _ 32 M^sMbhMt 



(A21) 



From our two angular momentum relations, we can de- 
rive a relation linking the mass loss rate to the change 
in binary separation, and determine the time-averaged 
mass loss rate. To do so, we will make a very simple 
assumption: that mass loss takes place at the location 
where we predict Roche-lobe overflow to occur. 

Solving Eq. ljA7|l for the binary separation and assum- 
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FIG. 16: Critical mass ratio for determining the stability of 
mass transfer from a relativistic polytropic NS companion as a 
function of the NS compactness, defined in terms of the poly- 
tropic index n (bottom), or the adiabatic index F = 1 + 1/n. 
Mass transfer is unstable to the right of the curve, stable to 
the left. The heavy solid lines depict the maximum possible 
compaction for a given polytropic EOS. The square and trian- 
gle represent the positions of the two models we are evolving 
dynamically. 



ing Rr = i?NS: we can insert Eq. ljA18|l and find 
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Plugging this result into Eq 
9-4n 



201, we find 
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(A23) 

Re-expressing the radiation-reaction angular momentum 
loss rate, Eq. I|A21(I . for a binary located at the Roche- 
lobe overfiow point, we find 
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We are now in a position to analyze the mass-transfer 
process from its onset until the final fate of the binary. 
To do so, we will make a few simplifying assumptions. 



First, we assume the NS is much less massive than the 
BH, such that Mbh ~ Mt and q <^ 1. Also, we will 
ignore the relativistic changes to the Newtonian power- 
law mass-radius relation, which has the effect of setting 
f — 1 uniformly (and thus df /dC = 0). In general, rel- 
ativistic polytropes expand more than their Newtonian 
analogues during mass loss, which increases the amount 
of orbital expansion seen for a given amount of mass loss. 
This statement would imply that a given amount of mass 
loss leads to a greater loss of angular momentum, or con- 
versely, that for a fixed angular momentum loss rate we 
would see a slight suppression of the mass loss rate. Un- 
der these assumptions, we find 
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where we define tcwo as the infall timescale when the 
NS first reaches the onset of mass transfer, as defined by 
Eq. Q, and where A/ns(0) is the NS mass at this time. 
This equation has the solution 
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Considering some familiar polytropic EOS, we find: 



n = 1.5 (r= 5/3) 
n = 1.0 (r = 2) 
n = 0.5 (r = 3) 
n = (r cx)) 
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Eq. (jA29p recovers the special value found by [43 , in the 
limit q <^ 1. We note that the typical rate found here 
satisfies A^ns ~ M^^/tawo [cf., Eq. (jA25p ]. while our 
dynamical simulations for cases in which viscosity is not 
important show Mns 3> Mns/towo- 



APPENDIX B: SYMMETRIES 

When constructing quasi-equilibrium configurations in 
either Newtonian or relativistic gravity, it is important 
to take note of the various symmetries present in the 
relevant equations. These symmetries can either be en- 
forced numerically, to save computational resources, or 
be used as a check to make sure that a numerical code is 
producing physically valid results. 

For quasi-equilibrium binary BHNS systems, virtually 
all gravitational formalisms will produce a configuration 
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that is equatorially symmetric so long as the NS spin 
axes are parallel to the orbital angular momentum axis, 
regardless of whether the NS is corotating, counterrotat- 
ing, or irrotational. If we fix the z-axis to be parallel to 
the various angular momenta mentioned above, we find 
that the transformation z —z, with a corresponding 
reflection for all vector and tensor quantities, leaves the 
hydrostatic equations and the metric invariant, and is 
compatible with the velocity field of the initial configu- 
ration as well. For quasi-equilibrium configurations eval- 
uated using full GR in the Kerr-Schild metric, this is the 
only symmetry plane of note present in the initial data. 
In our calculations, this symmetry is enforced at the code 
level: each SPH particle is treated as if it were a pair of 
particles, each of half the total mass, with one copy lying 
above the equatorial plane and one below. All spectral 
decompositions include this symmetry as well, setting to 
zero all spherical harmonics that are incompatible with 
a fully symmetric description. 

We show here that there is an additional symme- 
try plane for the case of equilibrium BHNS binaries in 
conformally fiat gravity, because the formalism is time- 
symmetric. Directions are defined such that the axis 
of separation between the BH and NS lies along the x- 
axis, and the orbital angular momentum points in the 
z-direction. Thus, the y-direction represents the direc- 
tion of motion for each object. Reversing the time di- 
rection requires us to perform two operations in order 
to maintain invariance. First, we must invert all vector 
quantities, most notably the velocity and the shift vector. 
Second, we must perform an inversion in the y-direction, 
to account for the inversion of our angular coordinate 
(f) —(f). We find the following relations are compatible 
with the initial velocity field, and leave our hydrostatic 
equations invariant: 

Scalars : f{x,y,z,t) — f{x,—y,z,—t), (Bl) 
Vectors : [w^, u", u^](a;, z, i) = 

-[v-,-vy,v']{x,-y,z,-t). (B2) 

Evaluating the expressions above at t = yields the sym- 
metries in the y-direction for our initial data. This sym- 
metry can be extended to tensor quantities as well: two 
index tensor elements satisfy the relation Tij{x, —y, z) — 
(— l)^+^«ry(x, 2), where Ny is the number of "y" in- 
dices present for a given element. The additional fac- 
tor of —1 is necessary to describe the inversion prop- 
erties under time symmetry. Equatorial symmetry can 
be described in the same language, Tyj, {x,y,—z) = 
(-ir^T,,fe...(a;,2/,z). 

It is straightforward to check that these symmetries are 
compatible with all equations governing the construction 
of quasi-equilibrium NS binaries in CF gravity. To do 
so, we note that any tensorial operation on fields satis- 
fying these symmetries will maintain these symmetries, 
including gradients and inner products. 

First, we note that the matter sources in the equations 
for the conformal factor and lapse function, Eqs. (|25|l and 



(|26|l . E, P, and S, as scalars, are symmetric in both y 
and z. The same pattern follows for scalar densities and 
quantities like u'^. As an immediate consequence, the 
lapse and conformal factor share the same symmetries, 
since the only other source terms involve KijK^^ , itself a 
scalar. 

The equilibrium velocity field is v — Q x r = 
[—ily, fix, 0], which satisfies 

[v-,vy,v']{x, ~y, z) = -[v\ -yy, v^']{x, y, z). (B3) 

Equatorial symmetry is satisfied trivially. In the CF case, 
the BH component of the shift vector is zero, and also 
satisfies the symmetry relations trivially. Since the shift 
equation, Eq. H24|l depends only on the velocity field and 
other tensors, we can extend the symmetry to describe 
all quantities present in the equation. This statement 
cannot be made for the Kerr-Schild metric, and serves 
as the simplest example of how time-symmetry fails. For 
the Kerr-Schild case, 

/^BH-KS oc a;* — > 

W.py, I3']{x, -y, z) = r , -f3\ y, z).(B4) 

Thus, the BH contribution to the shift satisfies a differ- 
ent symmetry relation than the initial velocity field, and 
there is no global symmetry present. The physical rea- 
sons underlying this fact are discussed in more detail in 
Appendix E of [lO^ . 

APPENDIX C: PRESERVING STATIONARY 
EQUILIBRIUM DURING CF EVOLUTION 

We show here that the CF evolution equations main- 
tain strict stationary equilibrium for initial data satisfy- 
ing both the thin-sandwich equations and the integrated 
Euler equation. That is, if we have an initially uniformly 
rotating matter configuration that satisfies the condition 
h/u^ = C, where C is constant throughout space at 
t = 0, and if the gravitational field satisfies the CTS 
equations, Eqs. (|24|l . H25I) . and (|26|l . the time derivatives 
of all quantities go to zero. Note that under these as- 
sumptions, this statement applies to both Eulerian and 
Lagrangian time derivatives, related by the expression 
d/dt = d/dt + v^d/dx^ since the difference term goes to 
zero when — 0, which applies in the corotating frame. 

First, we note that the continuity equation, Eq. I|15|l. 
and energy equation, Eq. H22I) . are trivially conserved 
when = 0, implying that and are conserved au- 
tomatically. 

The only evolution equation that requires a more thor- 
ough look is the Euler equation, Eq. (|19|l . which we re- 
express for convenience in the equivalent form 

We can show that the RHS of this expression is zero 
by starting from the integrated Euler equation, which 
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implies that 



0, 



(C2) 



From the relativistic Gibbs-Duhem relation, we know 
that 



d,h d,P 



(C3) 



The gradient of u° is determined from the normaUzation 
UaU"' = uqvP 

= 0. From this, we conclude 



(C4) 



and differentiating. 



■ 

Inserting the expression Ui = g^ivP = —ij^'^vPSijP^ for a 
configuration with = 0, we see that 

= - ^ahu'^d.^a + uAP'^ - j^d,^ j , (C6) 

and combining terms, 

^'-f , f u On , - a afe 2{tfc{tfe 

= 7; + ahu'dia + UkOifJ - , „ .. diip 



dui 
'~dt 



= 0. 



(C7) 



Thus, the RHS of the Euler equation, Eq. H19() . is zero 
under these assumptions, and Ui is also invariant. 



We now consider the field equations for the GTS initial 
data. The fields are described by a set of linked elliptic 
equations, Eqs. (^5)) . and whose source terms 

involve the fields themselves, as well as three quantities: 
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(G8) 

(G9) 
(CIO) 



where po, e, and P are the standard relativistic mass den- 
sity, internal energy density, and pressure, respectively. 
The Lorentz factor 7„ = au^, can be solved implicitly, 
from Eq. (|18|l . Goupled with our field equations, we have 
a set of six completely linked equations for six variables 
(7„, a, ip, (3^) and our conserved matter quantities (p*, 
e*, Ui). So long as a unique solution exists for our choice 
of matter variables, we know that this solution will re- 
main invariant so long as we choose our elliptic equation 
boundary conditions to be invariant. Thus, the RHS of 
all our evolution equations will remain zero, and the mat- 
ter configuration will remain in equilibrium. 
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